计算方法上机报告174四阶龙格-库塔法求解常微分方程的初值问题4.1算法原理及程序框图一阶常微分方程(组)和高阶常微分方程的初值问题最终都可以转化为一阶常微分方程组的初值问题,其向量形式为式(17)。0,,xxxaxbayfyyy(17)式(17)在形式上与单个微分方程的初值问题完全相同,只是将数量函数变成了向量函数。因此,标准四级四阶龙格—库塔法的向量形式和分量形式分别为式(18)和式(19)。1123411223431226,,22,22,iiiiiiiiiihxhhxhhxhxhyyKKKKKfyKKfyKKfyKfyK(18),1,12341121112121221222312411322331226;,,,;,,,1,2,,2222;,,,2222;,,,jijijjjjjiiininjiiininjiiinijiiininyyKKKKKhfxyyyKKKhKhfxyyyjnKKKhKhfxyyyKhfxhyKyKyK(19)根据式(19)可以得到使用四阶龙格-库塔法求解常微分方程组的初值问题的程序框图如图19。4四阶龙格-库塔法求解常微分方程的初值问题18图19标准四级四阶龙格—库塔法程序框图4.2程序使用说明标准四级四阶龙格—库塔法求解常微分方程初值问题的matlab程序见附录5,该程序可以用来求解一阶常微分方程(组)和高阶常微分方程的初值问题,运行程序按照命令窗口的提示输入相关变量直至得到结果。运行该程序会出现提示“请输入一阶微分方程组的个数或者高阶方程组的阶数n:”,此时需要输入一阶常微分方程组的个数,若求解的问题为高阶常微分方程组的初值问题,则需要将该高阶常微分方程转化为一阶常微分方程组,并输入阶数n。按Enter进行下一步则出现提示“请输入求解区间下限:”,此时需要输入区间[a,b]中的a值,输入完成后继续按Enter出现提示“请输入求解区间上限:”,根据提示输入b的值。继续按Enter出现提示“请输入求解步长h:”,需要输入计算步长。继续按Enter会出现提示“请依次输入n个方程的右端函数fi(x,y(x)):”此时需要依次输入式(17)中向量函数f(x,y(x))的各个元素,其中y(x)的各个分量的输入格式为y(1),y(2)……,例如,含有两个方程的一阶常微分方程组的第二个方程为2212sin22xyexyy,则需要输入exp(2*x)*sin(x)-2*y(1)+2*y(2)。输入完成之后按Enter出现提示“请依次输入n个初值yi(0):”此时则需要依次输入式(17)中的向量y0中的各元素。至此,所有输入已完成,按Enter开始进行计算,最终的结果储存在y矩阵中,同时将会出现y1随x的变化图像。计算方法上机报告194.3算例及计算结果《数值分析》课本第304页的例题计算实习9.2,用标准4级4阶R-K法求解23(0)1,(0)3,(0)2yyyyxyyy取步长h=0.05,计算y(1)的近似值,并与解析解y(x)=xex+2x-1作比较。首先,将高阶常微分方程转化为一阶常微分方程组1223332112323(0)1,(0)3,(0)2yyyyyyyyxyyy使用标准四级四阶龙格—库塔法求解常微分方程初值问题的matlab程序对以上一阶常微分方程组进行求解,计算结果如图20和图21所示。图204阶R-K法求解9.2的运行结果4四阶龙格-库塔法求解常微分方程的初值问题20图21计算得到的x~y1图像4阶R-K法数值解与解析解的对比如图22所示,从图中可以看出,4阶R-K法得到的数值解与解析解吻合的相当好。图22数值解与解析解的比较00.20.40.60.81-1-0.500.511.522.533.5xy1