matlab解龙格库塔,matlab龙格-库塔方法解微分方程.docx

matlab龙格-库塔方法解微分方程.docx

龙格库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格库塔方法为例，计算格式如下112342132436,,NNNNHYKFXYHKFXKY1龙格库塔法解一阶ODE对于形如的一阶ODE初值问题，可以直接套用公式，如今可以0,DYFXAB借助计算机方便的进行计算，下面给出一个实例2101DYX取步长H01，此处由数学知识可得该方程的精确解为。在这里利用MATLAB12YX编程，计算数值解并与精确解相比，代码如下1写出微分方程，便于调用和修改FUNCTIONVALODEFUNX,YVALY2X/YEND2编写RUNGEKUTTA方法的函数代码FUNCTIONYRUNGE_KUTTAH,X0,Y0K1ODEFUNX0,Y0K2ODEFUNX0H/2,Y0H/2K1K3ODEFUNX0H/2,Y0H/2K2K4ODEFUNX0H,Y0HK3YY0HK12K22K3K4/6END3编写主函数解微分方程,并观察数值解与精确解的差异CLEARALLH01X00Y01X01H1Y1RUNGE_KUTTAH,X0,Y0FORK1LENGTHXXKX0KHYK1RUNGE_KUTTAH,XK,YKENDZSQRT12XPLOTX,Y,’’HOLDONPLOTX,Z, R 结果如下图，数值解与解析解高度一致2龙格库塔法解高阶ODE对于高阶ODE来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。5027502YY这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。初始条件为，将方程降阶，引入一个向量型变量Y0Y故有Y20750YDYYT记则至此，二阶方程降阶为一阶1Y150Y方程组。值得注意的是此时再用龙格库塔法进行求解时，代入的将是一个Y向量。同样利用MATLAB进行计算，步长H005，时间周期为0,201编写ODE函数FUNCTIONYODEFUN1,Y0此处Y0为一个列向量，因为时间T未显含在一阶方程组中所以ODE函数的第一个参数为空，要根据具体情况而定。YY022000200Y02750Y01/500END2编写RUNGEKUTTA函数FUNCTIONYRKFAH,T0,Y0K1ODEFUN1T0,Y0K2ODEFUN1T0H/2,Y0H/2K1K3ODEFUN1T0H/2,Y0H/2K2K4ODEFUN1T0H,Y0HK3YY0HK12K22K3K4/6END3编写主函数CLEARALLH005T005H20T00Y000初值YCELL1,LENGTHTY{1}RKFAH,T0,Y0ZZEROS2,LENGTHTFORK1LENGTHTY{K1}RKFAH,T0,Y{K}Z1,KY{K}1Z2,KY{K}2ENDPLOTT,Z1,, R 位移Y的图像HOLDONPLOTT,Z2,速度Y’的图像求解结果如下图