四阶龙格库塔法解微分方程
一、四阶龙格库塔法解一阶微分方程
①一阶微分方程：cos
，初始值y(0)=0，求
y t
解区间[0 10]。
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));
for i=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
②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。
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;
for i=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=[t0 tf];
[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);
subplot(212)
plot(x_ode45,y_ode45,'--');
xlabel('t');
ylabel('x');
legend('Exact');
图2
二、四阶龙格库塔法解二阶微分方程
①二阶微分方程：
cos y t ，初始值y(0)=0，y'(0)=-1，求解区间[0 10]。
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));
for i=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
②二阶微分方程：
7.35499*cos y x ，初始值y(1)=0，y'(1)=0，求解
区间[0 25]。
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;
for i=1:length(x)-1
K1=RK_z(i);
L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 and L1
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); % K4 and L4 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('Variable x');
ylabel('Variable y');
A=[x,RK_y];
[y,T]=max(RK_y);
legend('RK4方法');
fprintf('角度峰值等于%d',y) %角度的峰值也就是π
fprintf('\n')
fprintf('周期等于%d',T*h) %周期
function w=rightf_sys1(x,y,z)
w=7.35499*cos(y);
end
图4
注：四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。