《matlab龙格-库塔方法解微分方程》由会员分享，可在线阅读，更多相关《matlab龙格-库塔方法解微分方程（6页珍藏版）》请在金锄头文库上搜索。
1、龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下112342132436,+nnnnhyKfxyhKfxKy1 龙格-库塔法解一阶 ODE对于形如 的一阶 ODE 初值问题，可以直接套用公式，如今可以0, dyfxab借助计算机方便的进行计算，下面给出一个实例2 101dyx取步长 h=0.1，此处由数学知识可得该方程的精确解为 。在这里利用 MATLAB12yx编程，计算数值解并与精确解相比，代码如下：(1)写出微分方程，便于调用和修改function val = odefun( x,y
2、 )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 allh = 0.1;x0 = 0;y0 = 1;x = 0.1:h:1;y(1) = runge_
3、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);plot(x,y,*);hold onplot(x,z,r);结果如下图，数值解与解析解高度一致2 龙格-库塔法解高阶 ODE对于高阶 ODE 来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。5027502yy&这是一个二阶 ODE,描述的是一个物体的有阻尼振动情况。初始条件为 ，将方程降阶，引入一个向量型变量 Y0 y&故有 Y20750ydYyt&记 则 至
4、此，二阶方程降阶为一阶1 y&150Y方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个 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,
5、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;t0 = 0;Y0 = 0;0;%初值Y = cell(1,length(t);Y1 = 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);%位移y的图像 hold onplot(t,z(2,:);%速度y 的图像求解结果如下图
