如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法Word下载.docx
《如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法Word下载.docx》由会员分享，可在线阅读，更多相关《如何用matlab求解非线性微分方程组基于龙格库塔的数值微分算法Word下载.docx（21页珍藏版）》请在冰豆网上搜索。
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
g_f_N（1,2）
g_f=g_f_space（1,x）;
fory=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
画出来的图形结果如下：
选择x和y轴都是log显示的模式，而原论文的图为：
这种代码的问题是，可能由于L不够大（L太大计算机太久，还没有算过，所以用可能），导致结果再tau_q比较大的时候，出现了问题，和原论文对不上。
如果用matlab自带的ode45或者ode23命令，也可以算，但是精度和老外PRL论文里的图形和数据差太多。
使用ode45命令的代码如下（程序名为Fig2_a_prl_v2.m）：
omega0=1;
%
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;
tau2=tau_space（1,2）;
myfun2a_v2'
[0tau2],[10]'
Er（1,2）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau3=tau_space（1,3）;
myfun3a_v2'
[0tau3],[10]'
Er（1,3）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau4=tau_space（1,4）;
myfun4a_v2'
[0tau4],[10]'
Er（1,4）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau5=tau_space（1,5）;
myfun5a_v2'
[0tau5],[10]'
Er（1,5）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau6=tau_space（1,6）;
myfun6a_v2'
[0tau6],[10]'
Er（1,6）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau7=tau_space（1,7）;
myfun7a_v2'
[0tau7],[10]'
Er（1,7）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau8=tau_space（1,8）;
myfun8a_v2'
[0tau8],[10]'
Er（1,8）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau9=tau_space（1,9）;
myfun9a_v2'
[0tau9],[10]'
Er（1,9）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau10=tau_space（1,10）;
myfun10a_v2'
[0tau10],[10]'
Er（1,10）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau11=tau_space（1,11）;
myfun11a_v2'
[0tau11],[10]'
Er（1,11）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau12=tau_space（1,12）;
myfun12a_v2'
[0tau12],[10]'
Er（1,12）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau13=tau_space（1,13）;
myfun13a_v2'
[0tau13],[10]'
Er（1,13）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau14=tau_space（1,14）;
myfun14a_v2'
[0tau14],[10]'
Er（1,14）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau15=tau_space（1,15）;
myfun15a_v2'
[0tau15],[10]'
Er（1,15）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau16=tau_space（1,16）;
myfun16a_v2'
[0tau16],[10]'
Er（1,16）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau17=tau_space（1,17）;
myfun17a_v2'
[0tau17],[10]'
Er（1,17）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau18=tau_space（1,18）;
myfun18a_v2'
[0tau18],[10]'
Er（1,18）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau19=tau_space（1,19）;
myfun19a_v2'
[0tau19],[10]'
Er（1,19）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau20=tau_space（1,20）;
myfun20a_v2'
[0tau20],[10]'
Er（1,20）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau21=tau_space（1,21）;
myfun21a_v2'
[0tau21],[10]'
gf=gc-10^（-1）;
Er（1,21）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau22=tau_space（1,22）;
myfun22a_v2'
[0tau22],[10]'
Er（1,22）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau23=tau_space（1,23）;
myfun23a_v2'
[0tau23],[10]'
Er（1,23）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau24=tau_space（1,24）;
myfun24a_v2'
[0tau24],[10]'
Er（1,24）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
tau25=tau_space（1,25）;
myfun25a_v2'
[0tau25],[10]'
Er（1,25）=omega0*（abs（y（y_final,2）））^2-（omega0*gf^2）/4*（abs（y（y_final,1）+y（y_final,2）））^2-...
figure;
plot（tau_space,Er）;
子程序代码如下：
myfun1a_v2.m
functiondy=myfun1a_v2（t,y）
dy=zeros（2,1）;
%
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））;
myfun2a_v2.m
functiondy=myfun2a_v2（t,y）%
升级会员 