微分方程求解

1、微分方程求解一、实验目的与要求1掌握用Matlab求微分方程及其方程组解的方法；2学会求微分方程近似解的欧拉折线法；3学会建立一些简单问题的微分方程模型，并能运用Matlab分析研究这些 问题。二、问题描述 对于很多实际问题，要直接找出所需的函数关系往往非常困难，但根据实 际问题所提供的条件， 有时却可以列出含有未知函数导数的关系式， 这样的关系 式就是所谓的微分方程。 怎样利用微分方程求得所需未知函数， 往往是我们解决 实际问题经常需要面对的问题，即解微分方程。这里我们借用Matlab对此问题 进行简单探讨。三、问题分析 在处理关于微分方程的实际问题时，我们一般须先建立微分方程，再利用 所学

2、的数学知识解微分方程。 事实上真正能找到精确解的微分方程只是很少一部 分，大部分只能求近似解，即数值解。四、试验过程1 求微分方程解析解的命令。求微分方程解析解的命令为：dsolve（方程1方程2，初始条件1初始条件2，自变量）, 对于可用积分方法求解的微分方程和微分方程组，可以用dsolve命令来求其通解和特解。例1：要求方程y 3y 4y 0的通解，可以输入以下语句Matlab命令：dsolve (D2y+3*Dy-4*y=0, x)ans =C1*exp(-4*x)+C2*exp(x)运行结果：即y CieC2ex4x推。如果自变量没有选定，默认自变量为2求微分方程数值解。求微分方程数值

3、解命令为ode45，ode23, ode15&对于不可以用积分方法求解的微分方程初值问题，可以用ode45，ode23,ode15s命令求特解。例4:求微分方程y yx y x2, yx 03的近似解(0 x 4)可用下 面的命令：fun cti on f=odef un 1(x,y)f=y*x+y+xA2;注：求一阶用D表示，二阶导数用D2表示，三阶导数用D3表示，以此类例2:解方程y4y 5y 0运行结果:dsolve (D2y+4*Dy+5*y=0,x)ans =C1*exp(-2*x)*si n(x)+C2*exp(-2*x)*cos(x)即:y e2xC1sin x C2co

4、sx例3:解方程1x2y 2xy xe运行结果:即:dsolve (1+xA2)*Dy+2*x*y=x*exp(xA2),x)ans=(1/2*exp(xA2)+C1)/(1+xA2)1x2-exC1y 2y1 x2如果要求微分方程的初值问题：y 4y 2y 0 yx 06 , y输入以下语句dsolve (D2y+4*Dy-2*y=0,y(0)=6,Dy(0)=10,x)运行结果：ans =x010，可(3+11/6A(1/2)*exp(-2+6A(1/2)*x)+(-11/6A(1/2)+3)*exp(-(2+6A(1/2)*x)即:y 311w61126 x3 e6x,y=ode45(o

5、defu n1,0,4,3);plot ( x , y , r -)输出结果为：例5:解初值问题y y t 1, y 01dsolve(Dy=-y+t+1,y(0)=1)输出结果：ans =t+exp(-t)即:y t et现在我们用数值求解命令求解后和解析解比较function f=odef un 2(t,y)f=-y+t+1;t=0:0.1:1; y=t+exp(-t); plot(t,y,b-)hold ont,y=ode45(odefu n2,0,1,1);plot(t,y,r.)hold off输出结果:例6:求初值问题y y sin 2x 0 , y 1 , y 1解：设yiy ,

6、 y2y，则原方程可化为yiy2y2y1sin 2xyi1 , y21Matlab语言：fun cti on f=odef un 3(x,y)f=y(2)；-y(1)-si n(2*x);x,y=ode45(odefu n3,pi,2*pi,1,1);plot(x,y(:,1),r-)输出结果:利用ode45命令还可以求解耦合微分方程，所谓耦合微分方程，方程组中 的未知函数是相互影响的，相互依赖的，其中的一个求解会影响到另一个求解， 下面求一对耦合微分方程的数值解：例7:解方程组xtyt其中x 00 , y 02.1y t 0.01y t sin x tMatlab语句：function f=

7、odef un 4(t,y)f=y (2),-0.01*y (2)-si n( y(1);t,y=ode15s(odefu n4,0,100,0,2.1);函数x x t的图像：plot(t,y(:,1),r-)输出结果为:函数y y t的图像：plot(t,y(:,2),r-)输出结果：用xt , yt生成参数图形plot(y(:,1),y(:,2),r-)输出结果:2iQ510152D253欧拉折线法对于初值问题y f x, y , y X。y，我们考虑函数y x的线性近似Lx y Xoy x。x x。由于函数y x可微，在包含X。的一个很小的邻域内L x是y x得很好的近 似。欧拉折线法

8、就是通过一系列的线性近似得到在较大区间内的yx的近似解。第一步：设xix0X，其中X很小，贝UyiL xiyoy xoXixoy。fxo, y。Xixo是y Xi得很好的近似，在区间xo, xi（无妨设x O） 上y x能被L x很好的近似。第二步：利用Xi, yi和斜率f & , yi来进行下一步近似，设X2Xix，y X2由y2yif Xi, yiX2Xi近似表示。第三步：利用点x2, y2和斜率f x2, y2，对于x3x2x，yx3由y3y2f X2, y2X3X3近似表示。 接这个点列的折线就是初值问题y f x, y , y x0y0的一个近似解。这就是 所谓的欧拉折线法。

9、这样我们就得到一列点列x0, y0, x1, y1, x2, y2, x3, y3而连其一般的步骤是x1x0dx ,y1y0f x0, y0 x1x0 x2x1dx ,y2y1f x1, y1x2x1x3x2dx ,y3y2f x2, y2x3x2xnxn1dx, ynyn1f xn 1, yn1xnxn 1例8： 利用欧拉折线法球初值问题y 1 y , y 0 1的近似解。以下是求此初值问题的Matlab语句function odefun6(n,d)X=0,1;for k=1:n/dX(k+1,1)=X(k,1)+d;X(k+1,2)=X(k,2)+(1+X(k,2)*d;endplot(X

10、(:,1),X(:,2)odefun6(1,0.01)此初值问题的精确解为y 2ex1，以上语句可以实现对精确解和近似解 的图像进行比较。hold onx=0:0.01:1;y=2*exp(x)-1;plot(x,y,r-)hold off输出结果:五、结论与应用研究卫星绕地球运行的轨迹。根据牛顿第二运动定律：F mam吐和万有引力定理Fdt2G啤。所以ra G鸟，其中M为地球的质量，x, yr为卫星所在位置的坐标，r . x2y2。因此我们有在x轴上加速度分量为aGx，r在y轴上加速度分量为ayGy，设卫星的运动方程为rxtyt假定卫星以初速度vy04000m/s在x 0,则有MG飞xr。如

11、果我们MGryr4.2 107m处入轨，地球质量为24M 5.97 10 kgG 6.672 1011N m2/kg2，这是一个初值问题。设y1x , y2y , yax , y4y微分方程可化为:yiy3Matlab语句:functionf=odefun5(t,y,flag,G,M)r=sqrt(y(1)A2+y(2)A2);f=y(3),y(4),-G*(M/rA3)*y(1),-G*(M/rA3)*y(2);G=6.672e-11;M=5.97e24;t,y=ode45(odefun5,0,60*60*24*6,-4.2e7,0,0,4000,G,M);plot(y(:,1),y(:,2

12、),r-)hold onX,Y,Z=sphere(10);axis(image)R=0.64e7;X=R*X;Y=R*Y;z=O*Z;surf(X,Y,Z,FaceColor,red,EdgeColor, non e);camlightright;phonghold off输出结果:y2y4y3y4yi0GMrGMryiy24.2 107, y200, yi0 0, y204000lighti ngKm六、练习求下列微分方程的通解。2)y 2y5yexsin 2x3)y 6y9y/ 3xx 1 e2.求初值问题yy sin2x 0 , y3.求微分方程x21 y 2xy cosx41)2y0yyx拉折线法得到的近似解0 x1 ,y1的解。0在初始条件yx。1下的精确解和用欧1并作图。注意观察dx的选取对解的精确度的影响。4. 求微分方