四阶龙格库塔法解微分方程

一、四阶龙格库塔法解一阶微分方程

①一阶微分方程： ，初始值y(0)=0，求解区间[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

②一阶微分方程： ，初始值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

二、四阶龙格库塔法解二阶微分方程

①二阶微分方程： ，初始值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

②二阶微分方程： ，初始值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

注：四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。