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

1、龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展 开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例, 计算格式如下y n 1KiK2K3K4VnKi2K2 2K3K4f Xn, Vnh hf XnVnh 2K1Xn22Kh - 2+ynf Xnh* hK1龙格-库塔法解一阶ODEdy f x,y a x对于形如dxb的一阶ODE初值问题，可以直接套用公式，如今可以y a V。借助计算机方便的进行计算，下面给出一个实例dydx取步长2xy 0 x 1 y1h=0.1 ,此处由数学知识可得该方程的精确解为y 小2x 。在这里利用 MATLAB编程，计算数

2、值解并与精确解相比，代码如下:(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)编写主函数解微分方程,并观察数值解与精确解的差异

3、clear allh = 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);endz = sqrt(1+2*x););piot(x,y hold onplot(x,z, r);结果如下图，数值解与解析解高度一致 Fig Lire 1 X文件内编事(E)查看(V)播入6 TMCB桌面葡口叫砰助(H)口直口通|0国白要，| 臣1|IQ2龙格-库塔法解高阶ODE对于高阶ODE来说，通用的方法是将高阶方程通过引

4、入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。500& 200/ 750 y 2000这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。初始条件为y 0 & 0 ,将方程降阶，引入一个向量型变量Yy&y 、& dY&Y。故有浴工T82000 200/ 750y&dty& -500Y 2记y Y 1 & Y 2则,2000 200Y 2750Y 1 至此，二阶方程降阶为一阶500方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个Y向量。同样利用MATLAB进行计算，步长 h=0.05,时间周期为0,20.(1)编写ODE函数function Y = odefun1(

5、,Y0 )%此处Y0为一个列向量，因为时间t未显含在一阶方程组中%所以。de函数的第一个参数为空，要根据具体情况而定。Y = 丫0(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 allh = 0.05;t = 0.05:h:20;Y 0 = 0;Y 0 = 0;0;%初值Y = cell(1,length(t);Y 1 = rkfa( h,t0,Y0 );z = zeros(2,length(t);for k=1:length(t)Yk+1 = rkfa( h,t0,Yk);z(1,k) = Yk(1);z(2,k) = Yk(2);endplot(t,z(1,:),r);%B 移 y