CWG-微分方程求解.doc

微分方程求解一、 实验目的与要求1 掌握用Matlab求微分方程及其方程组解的方法；2 学会求微分方程近似解的欧拉折线法；3 学会建立一些简单问题的微分方程模型，并能运用Matlab分析研究这些问题。二、 问题描述对于很多实际问题，要直接找出所需的函数关系往往非常困难，但根据实际问题所提供的条件，有时却可以列出含有未知函数导数的关系式，这样的关系式就是所谓的微分方程。怎样利用微分方程求得所需未知函数，往往是我们解决实际问题经常需要面对的问题，即解微分方程。这里我们借用Matlab对此问题进行简单探讨。三、 问题分析在处理关于微分方程的实际问题时，我们一般须先建立微分方程，再利用所学的数学知识解微分方程。事实上真正能找到精确解的微分方程只是很少一部分，大部分只能求近似解，即数值解。四、 背景知识介绍1 求微分方程解析解的命令。求微分方程解析解的命令为：dsolve（方程1，方程2，初始条件1，初始条件2，自变量），对于可用积分方法求解的微分方程和微分方程组，可以用dsolve命令来求其通解和特解。例1：要求方程的通解，可以输入以下语句Matlab命令： dsolve (D2y+3*Dy-4*y=0,x)运行结果： ans =C1*exp(-4*x)+C2*exp(x)即 注：求一阶用D表示，二阶导数用D2表示，三阶导数用D3表示，以此类推。如果自变量没有选定，默认自变量为t。例2：解方程 dsolve (D2y+4*Dy+5*y=0,x)运行结果： ans =C1*exp(-2*x)*sin(x)+C2*exp(-2*x)*cos(x)即： 例3：解方程 dsolve (1+x2)*Dy+2*x*y=x*exp(x2),x)运行结果：ans =(1/2*exp(x2)+C1)/(1+x2)即： 如果要求微分方程的初值问题：，可输入以下语句 dsolve (D2y+4*Dy-2*y=0,y(0)=6,Dy(0)=10,x)运行结果： ans = (3+11/6(1/2)*exp(-2+6(1/2)*x)+(-11/6(1/2)+3)*exp(-(2+6(1/2)*x)即： 2求微分方程数值解。求微分方程数值解命令为ode45，ode23，ode15s。对于不可以用积分方法求解的微分方程初值问题，可以用ode45，ode23，ode15s命令求特解。例4：求微分方程的近似解（）可用下面的命令： function f=odefun1(x,y) f=y2*x+y+x2; x,y=ode45(odefun1,0,4,3); plot ( x , y , r - )输出结果为： 例5：解初值问题 dsolve(Dy=-y+t+1,y(0)=1)输出结果：ans =t+exp(-t)即： 现在我们用数值求解命令求解后和解析解比较function f=odefun2(t,y)f=-y+t+1;t=0:0.1:1; y=t+exp(-t); plot(t,y,b-)hold ont,y=ode45(odefun2,0,1,1);plot(t,y,r.) hold off输出结果：例6：求初值问题解：设，则原方程可化为Matlab语言：function f=odefun3(x,y)f=y(2);-y(1)-sin(2*x);x,y=ode45(odefun3,pi,2*pi,1,1);plot(x,y(:,1),r-)输出结果：利用ode45命令还可以求解耦合微分方程，所谓耦合微分方程，方程组中的未知函数是相互影响的，相互依赖的，其中的一个求解会影响到另一个求解，下面求一对耦合微分方程的数值解：例7：解方程组其中Matlab语句： function f=odefun4(t,y)f=y(2),-0.01*y(2)-sin(y(1);t,y=ode15s(odefun4,0,100,0,2.1);函数的图像：plot(t,y(:,1),r-)输出结果为：函数的图像：plot(t,y(:,2),r-)输出结果：用生成参数图形plot(y(:,1),y(:,2),r-)输出结果：五、 实验过程1欧拉折线法对于初值问题，我们考虑函数的线性近似 由于函数可微，在包含的一个很小的邻域内是得很好的近似。欧拉折线法就是通过一系列的线性近似得到在较大区间内的的近似解。第一步：设，其中很小，则 是得很好的近似，在区间（无妨设）上能被很好的近似。第二步：利用和斜率来进行下一步近似，设，由 近似表示。第三步：利用点和斜率，对于，由 近似表示。 这样我们就得到一列点列。而连接这个点列的折线就是初值问题的一个近似解。这就是所谓的欧拉折线法。其一般的步骤是 例8：利用欧拉折线法球初值问题的近似解。以下是求此初值问题的Matlab语句 function odefun6(n,d) X=0,1; for k=1:n/d X(k+1,1)=X(k,1)+d; X(k+1,2)=X(k,2)+(1+X(k,2)*d; end plot(X(:,1),X(:,2) odefun6(1,0.01)此初值问题的精确解为，以上语句可以实现对精确解和近似解的图像进行比较。 hold on x=0:0.01:1; y=2*exp(x)-1; plot(x,y,r-)hold off输出结果：2微分方程的斜率场例9：求一阶微分方的斜率场。一阶微分方程求斜率场的Matlab语句如下： function odefun7(inx,axx,iny,axy) a=(axx-inx)/0.1; b=(axy-iny)/0.1; z=inx,iny;m=1;s=1; for j=1:b for k=1:a m=m+1;z(k+s,1)=z(k,1)+0.1;z(k+s,2)=z(k,2)+(j-1)*0.1; end s=m; end for i=1:length(z) x=z(i,1);y=z(i,2); x1=x+0.1;y1=y+(y2*x+y)*0.1; plot(x,x1,y,y1,b-) plot(z(:,1),z(:,2),b.,Markersize,2) hold on end hold off odefun7(-2,2,-2,2)输出结果为六、 结论与应用研究卫星绕地球运行的轨迹。根据牛顿第二运动定律：和万有引力定理。所以，其中为地球的质量，为卫星所在位置的坐标，。因此我们有在轴上加速度分量为，在轴上加速度分量为，设卫星的运动方程为，则有。如果我们假定卫星以初速度在处入轨，地球质量为。，这是一个初值问题。设微分方程可化为：Matlab语句：function f=odefun5(t,y,flag,G,M)r=sqrt(y(1)2+y(2)2);f=y(3),y(4),-G*(M/r3)*y(1),-G*(M/r3)*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),r-)hold onX,Y,Z=sphere(10);axis(image)R=0.64e7;X=R*X;Y=R*Y;z=0*Z;surf(X,Y,Z,FaceColor,red,EdgeColor,none);camlight right; lighting phonghold off输出结果：七、 练习1 求下列微分方程的通解。1）2）3）2 求初值问题的解。3 求微分方程在初始条件下的精确解和用欧拉折线法得到的近似解并作图。注意观察的选取对