MATLAB龙格-库塔方法解微分方程

1龙格-库塔法解一阶ODE

对于形如 的一阶ODE初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例

取步长h=0.1，此处由数学知识可得该方程的精确解为 。在这里利用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来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。

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

初始条件为 ，将方程降阶，引入一个向量型变量Y

故有  

记则至此，二阶方程降阶为一阶方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个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’的图像

求解结果如下图