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

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

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

①一阶微分方程：

cos y

t & ，初始值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

②一阶微分方程：

()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]。