2018年常微分方程的数值解法.doc

常微分方程的数值解法 　　常微分方程的数值解法 　　????江南大学信计1203柯恒 　　一、前言 　　对于很多微分方程，我们很难求出解析解，这时我们需要采取数值手段求解。在数值分析这门课中，老师讲到????=??(??,??)的数值解法，我们采用了欧拉格式，后退欧拉格式，梯形公式，改进欧拉格式及4阶龙格-库塔方法求解。而老师没讲关于微分方程组和高阶微分方程的解法。现在我来粗虐的说一下微分方程组及高阶微分方程的数值解，这里主要是借助matlab中的相关函数进行计算并仿真出图像。 dy二、相关问题 　　初值问题 　　问题1，微分方程组 　　问题描述：给定一个微分方程，并且告诉我们初始状态。我们便可求出整个过程的解。 给定下面方程（L表示省略号）： 　　?dy?dx?f1(x,y1,y2,?,yn) 　　??dy2?f(x,y,y,?,y)?112ndx????????dyn?f(x,y,y,?,y)112n??dx 　　而且初始状态??′??(??0)已知记成y(0). 　　问题解决：四阶龙格库塔方法： 　　??????=????(????,??????,??????,……,??????); 　　?????? 　　?????? 　　??????????????=????(????+,??????+???????,??????+???????,……,??????+???????) ????????=????(????+,??????+???????,??????+???????,……,??????+???????) ????????=????(????+,??????+?????????????+???????,……,??????+???????), 　　????,??+????=????,??+?(??????+?????????+?????????+?????????) 通过龙格库塔方法，我们可以计算后面点的值，在matlab中调用ode45来实现四阶龙格库塔方法的调用。 　　应用举例： 　　有一个同步地球卫星，现在加速到4km/s进行变轨，试分析变轨后卫星运动的轨迹。已知地球的质量M=5.97e24,引力常量的值为6.672e-11 　　问题建模：我们已经知道行星的运动是在一个平面上的。以地球位置为原点，地球与卫星的连线为x轴建立坐标系。则由题目条件开始时速度x,y轴的分量分别是0,4000,位移在x，y轴的分量分别是4.2e7,0进入新轨道。通过牛顿第二定律及万有引力定律： 　　??=??????? 　　??=???????????设??1=??，??2=??，??3=?? ，??4=?? ，则方程可以化为线性方程组 　　y1’=y3; 　　y2’=y4; 　　y3’=-G*M*y1*(y1^2+y2^2)^(-3/2); 　　y4’=-G*M*y2*(y1^2+y2^2)^(-3/2); 　　下面使用matlab画出变轨后的图像： 　　%%%%%微分方程组的函数 　　function dy=f1(t,y) 　　dy=zeros(4,1); 　　G=6.672e-11; 　　M=5.97e24; 　　dy(1)=y(3); 　　dy(2)=y(4); 　　r=(y(1)^2+y(2)^2)^(3/2); 　　dy(3)=-G*M*y(1)/r; 　　dy(4)=-G*M*y(2)/r; 　　end 　　%%%%在工作窗口调用ode45函数，并画出图像 　　 [t,y]=ode45(‘f1’,[0,60*60*24*6],[4.2e7,0,0,4000]); 　　 x1=y(:,1)’; 　　 y1=y(:,2)’; 　　 plot(‘x1,y1,’r’); 　　 plot(x1,y1,’r’); 　　 title(‘卫星变轨后图像’); 　　 xlabel(‘x’); 　　 ylabel(‘y’); 　　 text(0,0,’地球位置’); 　　 hold on; 　　 x2=zeros(size(x1)); 　　 y2=zeros(size(y1)); 　　 plot(x2,y1); 　　 plot(x1,y2); 　　 hold off 　　在matlab中画出的图像如下所示： 　　　　通过上图知运动方程是一个椭圆。 　　问题2：高阶微分方程的数值解： 　　问题描述：给定一个n阶微分方程，且知道这个方程早某初始点处函数值和1~n-1阶导数值。求解这个方程。 　　y(n)(n?1)??f(x,y,y,?,y) 　　另y??y1,y???y2,?,y(n?1)?yn?1，则方程可以化成下面形式： 　　　　?dy?dx?y1 　　??dy1?y2?dx????? ?dyn?2?yn?1?dx?dy?n?1?f(x,y,y1,