微分方程求解

微分方程求解

一、 实验目的与要求

1．掌握用Matlab 求微分方程及其方程组解的方法； 2．学会求微分方程近似解的欧拉折线法；

3．学会建立一些简单问题的微分方程模型，并能运用Matlab 分析研究这些

问题。

二、 问题描述

对于很多实际问题，要直接找出所需的函数关系往往非常困难，但根据实际问题所提供的条件，有时却可以列出含有未知函数导数的关系式，这样的关系式就是所谓的微分方程。怎样利用微分方程求得所需未知函数，往往是我们解决实际问题经常需要面对的问题，即解微分方程。这里我们借用Matlab 对此问题进行简单探讨。 三、 问题分析

在处理关于微分方程的实际问题时，我们一般须先建立微分方程，再利用所学的数学知识解微分方程。事实上真正能找到精确解的微分方程只是很少一部分，大部分只能求近似解，即数值解。 四、 背景知识介绍

1． 求微分方程解析解的命令。 求微分方程解析解的命令为： dsolve （‘方程1’，‘方程2’，…，‘初始条件1’，‘初始条件2’，…，‘自变量’），

对于可用积分方法求解的微分方程和微分方程组，可以用dsolve 命令来求其通解和特解。

例1：要求方程043=-'+''y y y 的通解，可以输入以下语句Matlab 命令： dsolve (‘D2y+3*Dy-4*y=0’,’x ’)

运行结果：ans =C1*exp(-4*x)+C2*exp(x)

即x x e C e C y 241+=-

注：求一阶用D 表示，二阶导数用D2表示，三阶导数用D3表示，以此类推。如果自变量没有选定，默认自变量为’t ’。

例2：解方程054=+'+''y y y

dsolve ('D2y+4*Dy+5*y=0','x')

运行结果：ans =C1*exp(-2*x)*sin(x)+C2*exp(-2*x)*cos(x)

即： ()x C x C e y x c o s s i n 212+=- 例3：解方程(

)

2

212

x xe xy y x =+'+

dsolve ('(1+x^2)*Dy+2*x*y=x*exp(x^2)','x')

运行结果：ans =(1/2*exp(x^2)+C1)/(1+x^2)

即： 2

1

1212

x C e y x ++=

如果要求微分方程的初值问题：10,60240

='

==-'+''==x x y y

y y y ，可

输入以下语句

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)

即：()()x

62x 6261136113--+-???? ?

?-+???? ??+=e

e y 2．求微分方程数值解。

求微分方程数值解命令为ode45，ode23，ode15s 。

对于不可以用积分方法求解的微分方程初值问题，可以用ode45，ode23，ode15s 命令求特解。

例4：求微分方程3,0

22=++='=x y

x y x y y 的近似解（40≤≤x ）可用下

面的命令：

function f=odefun1(x,y) f=y^2*x+y+x^2; [x,y]=ode45(‘odefun1’,[0,4],3); plot ( x , y , ‘ r - ’ )

输出结果为：

例5：解初值问题()10,1=++-='y t y y dsolve('Dy=-y+t+1','y(0)=1')

输出结果：ans =t+exp(-t)

即：t e t y -+=

现在我们用数值求解命令求解后和解析解比较

function f=odefun2(t,y) f=-y+t+1;

t=0:0.1:1;y=t+exp(-t);plot(t,y,'b-') hold on

[t,y]=ode45('odefun2',[0,1],1); plot(t,y,'r.')

hold off

输出结果：

例6：求初值问题()()1,1,02sin ='==++''ππy y x y y 解：设y y y y '==21,，则原方程可化为

()()???

??==--='='1

,12sin 21

12

21ππy y x y y y y 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：解方程组

()()

()()()

()

?

?

?

-

-

=

'

=

'

t x

t y

t

y

t y

t

x

sin

01

.0

其中()()1.2

,

0=

=y

x

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]); 函数()t x

x=的图像：

plot(t,y(:,1),'r-')

输出结果为：

函数()t y

y=的图像：

plot(t,y(:,2),'r-')

输出结果：

用[][]t y t x ,生成参数图形

plot(y(:,1),y(:,2),'r-')

输出结果：

五、 实验过程 1．欧拉折线法

对于初值问题()()00,,y x y y x f y =='，我们考虑函数()x y 的线性近似

()()()()000x x x y x y x L -'+=

由于函数()x y 可微，在包含的一个很小的邻域内()x L 是()x y 得很好的近似。欧拉折线法就是通过一系列的线性近似得到在较大区间内的()x y 的近似解。

第一步：设x x x ?+=01，其中很小，则

()()()()()010********,x x y x f y x x x y y x L y -+=-'+==

是()1x y 得很好的近似，在区间[]10,x x （无妨设0>?x ）上()x y 能被()x L 很好的近似。

第二步：利用()11,y x 和斜率()11,y x f 来进行下一步近似，设x x x ?+=12，

()2x y 由

()()121112,x x y x f y y -+= 近似表示。

第三步：利用点()22,y x 和斜率()22,y x f ，对于x x x ?+=23，()3x y 由

()()332223,x x y x f y y -+= 近似表示。

………………

这样我们就得到一列点列()()()() 33221100,,,,,,,y x y x y x y x 。而连接这个点列的折线就是初值问题()()00,,y x y y x f y =='的一个近似解。这就是所谓的欧拉折线法。

其一般的步骤是

()()01000101,,x x y x f y y dx x x -+=+=()()

12111212,,x x y x f y y dx x x -+=+=()()

23222323,,x x y x f y y dx x x -+=+=………………

()()11111,,------+=+=n n n n n n n n x x y x f y y dx x x

例8：利用欧拉折线法球初值问题()10,1=+='y y y 的近似解。 以下是求此初值问题的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)

此初值问题的精确解为12-=x e y ，以上语句可以实现对精确解和近似解

的图像进行比较。

hold on

x=0:0.01:1;

y=2*exp(x)-1;

plot(x,y,'r-')

hold off

输出结果：

2．微分方程的斜率场

'2的斜率场。

例9：求一阶微分方y

=

y+

y

x

一阶微分方程求斜率场的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+(y^2*x+y)*0.1;

plot([x,x1],[y,y1],'b-')

plot(z(:,1),z(:,2),'b.','Markersize',2)

hold on

end

hold offodefun7(-2,2,-2,2)

输出结果为

六、 结论与应用

研究卫星绕地球运行的轨迹。

根据牛顿第二运动定律：22dt r

d m ma F ==和万有引力定理2r mM G F -=。所以

2

r

M G

a -=，其中为地球的质量，()y x ,为卫星所在位置的坐标，2

2y x r +=。因此我们有在轴上加速度分量为x r M

G a x 3-=，在轴上加速度分量为

y r M G a y 3-=，设卫星的运动方程为()()???==t y y t x x ，则有()()???

????

-=''-=''y r M G t y x r M G t x 33。如果我们

假定卫星以初速度()s m v y /40000=在()m x 7102.40?-=处入轨，地球质量为

kg M 241097.5?=。()

2211kg /m N 10672.6G ??=-，这是一个初值问题。设y y x y y y x y '='===4321,,,微分方程可化为：

()()()()????

?

?????

???='='=?-=-='-='='='40000,00,00,102.4021

27

123

413

34231y y y y y r

M G y y r M G y y y y y Matlab 语句：

function f=odefun5(t,y,flag,G ,M) r=sqrt(y(1)^2+y(2)^2);

f=[y(3),y(4),-G*(M/r^3)*y(1),-G*(M/r^3)*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 on

[X,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 phong hold off 输出结果：

七、 练习

1．求下列微分方程的通解。 1）()024=''+'''-y y y

2）x e y y y x 2sin 52=+'-'' 3）()x e x y y y 3196+=+'-'' 2．求初值问题1,1,02sin ='

==++''==π

π

x x y y

x y y 的解。

3．求微分方程()

0cos 212=-+'-x xy y x 在初始条件10

==x y

下的精确解和用欧

拉折线法得到的近似解()10≤≤x 并作图。注意观察的选取对解的精确度的影响。

4．求微分方程组???????==+-==-+==4

,081,0300t t y y x dt

dy x y x dt

dx

的特解，并作出用[][]t y t x ,生成参数

图形。

5．请用欧拉折线法在Matlab 下编写求解一阶微分方程初值问题的程序。

6．设河边点的正对岸点为：河宽h OA =，两岸为平行直线，水流速度为米每分有一只鸭子从点游向点，设鸭子在静水中的游动速度为米每分，且鸭子游动方向始终朝着点，求鸭子游过的迹线方程。并用Matlab 作出轨迹图。