四阶龙格库塔法解微分方程.docx
《四阶龙格库塔法解微分方程.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔法解微分方程.docx（8页珍藏版）》请在冰豆网上搜索。
四阶龙格库塔法解微分方程
四阶龙格库塔法解微分方程
一、四阶龙格库塔法解一阶微分方程
①一阶微分方程：
，初始值y（0）=0，求解区间[010]。
MATLAB程序：
%%%%%%%%%%%四阶龙哥库塔法解一阶微分方程
%%%%%%%%%%%y'=cost
%%%%%%%%%%%y（0）=0,0≤t≤10，h=0.01
%%%%%%%%%%%y=sint
h=0.01;
hf=10;
t=0:
h:
hf;
y=zeros（1,length（t））;
y
（1）=0;
F=@（t,y）（cos（t））;
fori=1:
（length（t）-1）
k1=F（t（i）,y（i））;
k2=F（t（i）+h/2,y（i）+k1*h/2）;
k3=F（t（i）+h/2,y（i）+k2*h/2）;
k4=F（t（i）+h,y（i）+k3*h）;
y（i+1）=y（i）+1/6*（k1+2*k2+2*k3+k4）*h;
end
subplot（211）
plot（t,y,'-'）
xlabel（'t'）;
ylabel（'y'）;
title（'Approximation'）;
span=[0,10];
p=y
（1）;
[t1,y1]=ode45（F,span,p）;
subplot（212）
plot（t1,y1）
xlabel（'t'）;
ylabel（'y'）;
title（'Exact'）;
图1
②一阶微分方程：
，初始值x
（1）=2，求解区间[13]。
MATLAB程序：
%%%%%%%%%%%四阶龙哥库塔法解微分方程
%%%%%%%%%%%x'（t）=（t*x-x^2）/t^2
%%%%%%%%%%%x
（1）=2,1≤t≤3,h=1/128
%%%%%%%%%%%精确解：
x（t）=t/（0.5+lnt）
h=1/128;%%%%%步长
tf=3;
t=1:
h:
tf;
x=zeros（1,length（t））;
x
（1）=2;%%%%%初始值
F_tx=@（t,x）（t.*x-x.^2）./t.^2;
fori=1:
（length（t）-1）
k_1=F_tx（t（i）,x（i））;
k_2=F_tx（t（i）+0.5*h,x（i）+0.5*h*k_1）;
k_3=F_tx（（t（i）+0.5*h）,（x（i）+0.5*h*k_2））;
k_4=F_tx（（t（i）+h）,（x（i）+k_3*h））;
x（i+1）=x（i）+（1/6）*（k_1+2*k_2+2*k_3+k_4）*h;
end
subplot（211）
plot（t,x,'-'）;
xlabel（'t'）;
ylabel（'x'）;
legend（'Approximation'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%ode45求精确解
t0=t
（1）;x0=x
（1）;
xspan=[t0tf];
[x_ode45,y_ode45]=ode45（F_tx,xspan,x0）;
subplot（212）
plot（x_ode45,y_ode45,'--'）;
xlabel（'t'）;
ylabel（'x'）;
legend（'Exact'）;
图2
二、四阶龙格库塔法解二阶微分方程
①二阶微分方程：
，初始值y（0）=0，y'（0）=-1，求解区间[010]。
MATLAB程序：
%%%%%%%%%龙格库塔欧拉方法解二阶微分方程
%%%%%%%%%y''=cost
%%%%%%%%%y'（0）=-1,y（0）=0
%%%%%%%%%转换为一阶微分方程
h=0.01;
ht=10;
t=0:
h:
ht;
%%%%%%%%%y'=z1,z1'=y''=cost
y=zeros（1,length（t））;
z=zeros（1,length（t））;
y
（1）=0;
z
（1）=-1;
f=@（t,y,z）（z）;
g=@（t,y,z）（sin（t））;
fori=1:
（length（t）-1）
K1=f（t（i）,y（i）,z（i））;
L1=g（t（i）,y（i）,z（i））;
K2=f（t（i）+h/2,y（i）+h/2*K1,z（i）+h/2*L1）;
L2=g（t（i）+h/2,y（i）+h/2*K1,z（i）+h/2*L1）;
K3=f（t（i）+h/2,y（i）+h/2*K2,z（i）+h/2*L2）;
L3=g（t（i）+h/2,y（i）+h/2*K2,z（i）+h/2*L2）;
K4=f（t（i）+h,y（i）+h*K3,z（i）+h*L3）;
L4=g（t（i）+h,y（i）+h*K3,z（i）+h*L3）;
y（i+1）=y（i）+h/6*（K1+2*K2+2*K3+K4）;
z（i+1）=z（i）+h/6*（L1+2*L2+2*L3+L4）;
end
plot（t,y）
xlabel（'t'）;
ylabel（'y'）;
title（'y''=sin（t）'）;
legend（'y''=sint'）;
图3
②二阶微分方程：
，初始值y
（1）=0，y'
（1）=0，求解区间[025]。
MATLAB程序：
%%%%%%%%%龙格库塔欧拉方法解二阶微分方程
%%%%%%%%%y''=7.35499*cosx
%set（0,'RecursionLimit',500）
h=0.01;
a=0;
b=25;
x=a:
h:
b;
RK_y
（1）=0;%初值
RK_z
（1）=0;
fori=1:
length（x）-1
K1=RK_z（i）;
L1=rightf_sys1（x（i）,RK_y（i）,RK_z（i））;%K1andL1
K2=RK_z（i）+0.5*h*L1;
L2=rightf_sys1（x（i）+0.5*h,RK_y（i）+0.5*h*K1,RK_z（i）+0.5*h*L1）;
K3=RK_z（i）+0.5*h*L2;
L3=rightf_sys1（x（i）+0.5*h,RK_y（i）+0.5*h*K2,RK_z（i）+0.5*h*L2）;
K4=RK_z（i）+h*L3;
L4=rightf_sys1（x（i）+h,RK_y（i）+h*K3,RK_z（i）+h*L3）;%K4andL4
RK_y（i+1）=RK_y（i）+1/6*h*（K1+2*K2+2*K3+K4）;
RK_z（i+1）=RK_z（i）+1/6*h*（L1+2*L2+2*L3+L4）;
end
plot（x,RK_y,'b-'）;
xlabel（'Variablex'）;
ylabel（'Variabley'）;
A=[x,RK_y];
[y,T]=max（RK_y）;
legend（'RK4方法'）;
fprintf（'角度峰值等于%d',y）%角度的峰值也就是π
fprintf（'\n'）
fprintf（'周期等于%d',T*h）%周期
functionw=rightf_sys1（x,y,z）
w=7.35499*cos（y）;
end
图4
注：
四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。
升级会员 