如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法.docx
- 文档编号:7295438
- 上传时间:2023-01-22
- 格式:DOCX
- 页数:21
- 大小:499.74KB
如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法.docx
《如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法.docx》由会员分享,可在线阅读,更多相关《如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法.docx(21页珍藏版)》请在冰豆网上搜索。
如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法
如何用matlab求解非线性微分方程组(基于龙格库塔的数值微分算法)?
例如要求解下面这个含时间的线性微分方程组,如下图所示
其中:
这里的tau_q,g_f,w0都是不含时间的常数:
初始条件为:
归一化条件:
通过观察以上方程组,可以看到:
每个方程等式里面并没有u(t)*v(t)这样的交叉项形式,也就是没有非线性项。
因此,用龙格库塔的数值微分算法,可以有两种思路,一种思路是利用哈密顿量矩阵去算,另外一种思路是用列向量去算。
两种代码如下解析:
(1)利用哈密顿量矩阵的思路
先构建一个子程序m文件,文件名为:
coupled_differential_equation.m
这里.m文件写入以下代码:
function[U,Er,U_condition,TEST]=coupled_differential_equation(rho_0,T,L,w0,g_f,tau_q)
t=linspace(0,T,L);%最终演化时间T就等于每一个tau_q
h=T/L;%h这是步长和精度,不同的tau_q,精度应该不一样
C{1}=rho_0;%演化矩阵的初态
forx=1:
1:
L
TEST(1:
2,x)=C{1};%这里的测试代码时用,没有其他用途,以下才是龙格库塔算法
k1=master(t(1,x),C{1},w0,g_f,tau_q);
k2=master(t(1,x)+h/2,C{1}+h*k1/2,w0,g_f,tau_q);
k3=master(t(1,x)+h/2,C{1}+h*k2/2,w0,g_f,tau_q);
k4=master(t(1,x)+h,C{1}+h*k3,w0,g_f,tau_q);
C{1}=C{1}+(k1+2*k2+2*k3+k4)*h/6;
U=C{1};
end
Er=w0*(abs(U(2,1)))^2-(w0*g_f^2)/4*(abs(U(1,1)+U(2,1)))^2-(w0*sqrt(1-g_f^2)-w0)/2;%这是根据Er公式算写的,U(1,1)和U(2,1)就是演化到每一个最终时间T=tau_q时的解u(tau_q)和v(tau_q)解。
U_condition=(abs(U(1,1)))^2-(abs(U(2,1)))^2;%用来验证这两个值是否归一化的
functiony=master(t,P,w0,g_f,tau_q)
H=zeros(2,2);
g_t=g_f/tau_q*t;
H=[w0/i*(1-1/2*g_t^2),-w0/i*1/2*g_t^2;
w0/i*1/2*g_t^2,-w0/i*(1-1/2*g_t^2);];%微分方程组对应第一列u(t),第二列v(t)的形式去改写,之后把系数矩阵写出来就是这个形式
y=H*P;
主程序Fig2_a_prl.m代码如下:
clear;
tic
w0=1;
rho_0=[1;0];
g_c=1;
L=200000;%L的大小(精度)随着tau_q的增加而迅速增加!
注意tau_q变化需要调试这里
tau_x1=-1;
tau_x2=5;
tau_N=25;
tau_space=logspace(tau_x1,tau_x2,tau_N);
test_Er=zeros(1,tau_N);
g_f_space=[g_c,g_c-10^(-4),g_c-10^(-3),g_c-10^(-2),g_c-10^(-1)];
g_f_N=size(g_f_space);
figure
forx=1:
1:
g_f_N(1,2)
g_f=g_f_space(1,x);
fory=1:
1:
tau_N
tau_q=tau_space(1,y);
T=tau_q;
[U,Er,U_condition,TEST]=coupled_differential_equation(rho_0,T,L,w0,g_f,tau_q);
test_Er(x,y)=Er;
end
holdon
plot(tau_space,test_Er(x,:
));
x
end
画出来的图形结果如下:
选择x和y轴都是log显示的模式,而原论文的图为:
这种代码的问题是,可能由于L不够大(L太大计算机太久,还没有算过,所以用可能),导致结果再tau_q比较大的时候,出现了问题,和原论文对不上。
如果用matlab自带的ode45或者ode23命令,也可以算,但是精度和老外PRL论文里的图形和数据差太多。
使用ode45命令的代码如下(程序名为Fig2_a_prl_v2.m):
clear;
omega0=1;
tau_x1=-1;
tau_x2=5;
tau_N=25;
tau_space=logspace(tau_x1,tau_x2,tau_N);
%
t=0;y=0;y_size=0;
tau1=tau_space(1,1);
[t,y]=ode45('myfun1a_v2',[0tau1],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,1)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau2=tau_space(1,2);
[t,y]=ode45('myfun2a_v2',[0tau2],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,2)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
t=0;y=0;y_size=0;
tau3=tau_space(1,3);
[t,y]=ode45('myfun3a_v2',[0tau3],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,3)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau4=tau_space(1,4);
[t,y]=ode45('myfun4a_v2',[0tau4],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,4)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau5=tau_space(1,5);
[t,y]=ode45('myfun5a_v2',[0tau5],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,5)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau6=tau_space(1,6);
[t,y]=ode45('myfun6a_v2',[0tau6],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,6)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau7=tau_space(1,7);
[t,y]=ode45('myfun7a_v2',[0tau7],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,7)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau8=tau_space(1,8);
[t,y]=ode45('myfun8a_v2',[0tau8],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,8)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau9=tau_space(1,9);
[t,y]=ode45('myfun9a_v2',[0tau9],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,9)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau10=tau_space(1,10);
[t,y]=ode45('myfun10a_v2',[0tau10],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,10)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau11=tau_space(1,11);
[t,y]=ode45('myfun11a_v2',[0tau11],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,11)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau12=tau_space(1,12);
[t,y]=ode45('myfun12a_v2',[0tau12],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,12)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau13=tau_space(1,13);
[t,y]=ode45('myfun13a_v2',[0tau13],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,13)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau14=tau_space(1,14);
[t,y]=ode45('myfun14a_v2',[0tau14],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,14)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau15=tau_space(1,15);
[t,y]=ode45('myfun15a_v2',[0tau15],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,15)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau16=tau_space(1,16);
[t,y]=ode45('myfun16a_v2',[0tau16],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,16)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau17=tau_space(1,17);
[t,y]=ode45('myfun17a_v2',[0tau17],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,17)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau18=tau_space(1,18);
[t,y]=ode45('myfun18a_v2',[0tau18],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,18)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau19=tau_space(1,19);
[t,y]=ode45('myfun19a_v2',[0tau19],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,19)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau20=tau_space(1,20);
[t,y]=ode45('myfun20a_v2',[0tau20],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc;
Omega=100*omega0;
Er(1,20)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau21=tau_space(1,21);
[t,y]=ode45('myfun21a_v2',[0tau21],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc-10^(-1);
Omega=100*omega0;
Er(1,21)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau22=tau_space(1,22);
[t,y]=ode45('myfun22a_v2',[0tau22],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc-10^(-1);
Omega=100*omega0;
Er(1,22)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau23=tau_space(1,23);
[t,y]=ode45('myfun23a_v2',[0tau23],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc-10^(-1);
Omega=100*omega0;
Er(1,23)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau24=tau_space(1,24);
[t,y]=ode45('myfun24a_v2',[0tau24],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc-10^(-1);
Omega=100*omega0;
Er(1,24)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
%
t=0;y=0;y_size=0;
tau25=tau_space(1,25);
[t,y]=ode45('myfun25a_v2',[0tau25],[10]');
y_size=size(y);
y_final=y_size(1,1);
gc=1;
gf=gc-10^(-1);
Omega=100*omega0;
Er(1,25)=omega0*(abs(y(y_final,2)))^2-(omega0*gf^2)/4*(abs(y(y_final,1)+y(y_final,2)))^2-...
(omega0*sqrt(1-gf^2)-omega0)/2;
figure;
plot(tau_space,Er);
子程序代码如下:
myfun1a_v2.m
functiondy=myfun1a_v2(t,y)
dy=zeros(2,1);
omega0=1;
gc=1;
gf=gc;%
Omega=100*omega0;
tau_x1=-1;
tau_x2=5;
tau_N=25;
tau_space=logspace(tau_x1,tau_x2,tau_N);
gt=gf*t/tau_space(1,1);
dy
(1)=omega0/i*((1-(gt^2)/2)*y
(1)-(gt^2)/2*y
(2));
dy
(2)=omega0/(-i)*((1-(gt^2)/2)*y
(2)-(gt^2)/2*y
(1));
end
myfun2a_v2.m
functiondy=myfun2a_v2(t,y)%
dy=zeros(2,1);
omega0=1;
gc=1;
gf=gc;
Omega=100*omega0;
tau_x1=-1;
tau_x2=5;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 何用 matlab 求解 非线性 微分 方程组 基于 龙格库塔 数值 算法