(完整word版)MATLAB龙格-库塔方法解微分方程

hyyK12K22K3K4n1n6K1fxn,ynhh K2fxn,ynK122hhK3fxn,yn+K222Kfxh,yhKnn34

1龙格-库塔法解一阶ODE

dyfx,y axb对于形如dx 的一阶ODE初值问题，可以直接套用公式，如今可以借助计

yay0算机方便的进行计算，下面给出一个实例

2xdyy 0x1dxy y01取步长h=0.1，此处由数学知识可得该方程的精确解为计算数值解并与精确解相比，代码如下：

y12x 。在这里利用MATLAB编程，

(1)写出微分方程，便于调用和修改 function val = odefun( x,y )

val = y-2*x/y; end

(2)编写runge-kutta方法的函数代码

function y = runge_kutta( h,x0,y0 ) k1 = odefun(x0,y0);

k2 = odefun(x0+h/2,y0+h/2*k1); k3 = odefun(x0+h/2,y0+h/2*k2); k4 = odefun(x0+h,y0+h*k3);

y = y0+h*(k1+2*k2+2*k3+k4)/6; end

(3)编写主函数解微分方程,并观察数值解与精确解的差异

clear all h = 0.1; x0 = 0; y0 = 1; x = 0.1:h:1;

y(1) = runge_kutta(h,x0,y0); for k=1:length(x) x(k) = x0+k*h;

y(k+1) = runge_kutta(h,x(k),y(k)); end

z = sqrt(1+2*x);

plot(x,y,’*’); hold on plot(x,z,'r');

结果如下图，数值解与解析解高度一致

2龙格-库塔法解高阶ODE

对于高阶ODE来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个

实例进行说明。

500y200y750y2000

这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。 初始条件为

y00 y00 ，将方程降阶，引入一个向量型变量Y

ydYyy2000200y750y Y故有Ydtyy500Y2记yY1 yY2则Y2000200Y2750Y1至此，二阶方程降阶为一阶方程

500组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个Y向量。同样利用MATLAB进行计算，步长h=0.05，时间周期为[0,20].

(1) 编写ODE函数

function Y = odefun1( ~,Y0 )

% 此处Y0为一个列向量，因为时间t未显含在一阶方程组中 % 所以ode函数的第一个参数为空，要根据具体情况而定。 Y = [Y0(2);

(2000-200*Y0(2)-750*Y0(1))/500;]; end

(2) 编写runge-kutta函数

function Y = rkfa( h,t0,Y0 ) k1 = odefun1(t0,Y0);

k2 = odefun1(t0+h/2,Y0+h/2*k1);

k3 = odefun1(t0+h/2,Y0+h/2*k2); k4 = odefun1(t0+h,Y0+h*k3); Y = Y0+h*(k1+2*k2+2*k3+k4)/6; end

(3) 编写主函数

clear all h = 0.05; t = 0.05:h:20; t0 = 0; Y0 = [0;

0];%初值 Y = cell(1,length(t)); Y{1} = rkfa( h,t0,Y0 ); z = zeros(2,length(t));

for k=1:length(t)

Y{k+1} = rkfa( h,t0,Y{k}); z(1,k) = Y{k}(1);

z(2,k) = Y{k}(2);

end

plot(t,z(1,:),'r');%位移y的图像

hold on

plot(t,z(2,:));%速度y’的图像

求解结果如下图