微分方程求解

微分方程的求解

----基于Maple工具一、Maple工具简介

Maple是加拿大滑铁卢大学（Waterloo University）研制的一种计算机代数系统。经过近20年的不断发展，数学软件Maple已成为当今世界上最优秀的几个数学软件之一，它以良好的使用环境、强有力的符号计算能力、高精度的数字计算、灵活的图形显示和高效的可编程功能,为越来越多的教师、学生和科研人员所喜爱，并成为他们进行数学处理的工具。可以容易的运用Maple软件解决微积分、解析几何、线性代数、微分方程、计算方法、概率统计等数学分支中的常见的计算问题。1980 年9月，加拿大Waterloo大学的符号计算研究小组成立，开始了符号计算在计算机上实现的研究项目。数学软件Maple是这个项目的产品。目前，这仍是一个正在研究的项目。Maple的第一个商业版本Maple3.3是1985年出版的。随着几经更新，Maple被广泛的使用，得到越来越多的用户。

Maple 软件主要由三部分组成：用户界面（Iris），代数运算器（kernel），外部函数库（External library）。用户界面和代数运算器是用C语言写的，只占整个软件的一小部分，当系统启动时，即被装入。Iris负责输入命令和算式的初步处理、显示结果、函数图像的显示等。Kernel负责输入的编译、基本的代数运算，如有理数运算、初等代数运算，还负责内存管理。Maple的大部分数学函数

和过程是用Maple自身的语言写成的，存于外部函数库中。当一个函数调用时，在多数情况下，Maple会自动将该函数的过程调入内存，一些不常用的函数才需要用户自己将它们调入。另外有一些特别的函数包也需要用户自己调入，如线性代数包、统计包，这使得Maple 在资源的利用上具有很大的优势，只有最有用的东西才留住内存，这是Maple可以在较小内存的计算机上正常运行的原因。

Maple不仅仅提供编程工具，更重要的是提供数学知识。用户在直观的文档计算环境中完成各种数学问题，无限精度的计算消除了计算误差。Maple是教授、研究员、科学家、工程师、学生们必备的科学计算工具，从简单的数字计算到高度复杂的非线性问题，Maple 都可以帮助您快速、高效地解决问题。用户可在单一的环境中完成多领域建模和模拟，符号计算、数值计算、程序设计、技术文件、报告演示、算法开发、与外部程序连接等功能，满足各个层次用户的需要，从高中学生到高级研究人员。

Maple 有三大关键特征：数学引擎、操作简单、外部连接性。

数学引擎：Mathematics＝Maplesoft ！解决数学问题时，世界上没有任何其他软件比Maple更完整、更好。

操作简单：Maple人性化的界面让用户只需要按几个键就可以解决大量复杂的计算问题，Maple的文件模式界面可以创建多样化的、专业级的技术文件，并可以自由转换为其他格式的文件，如Latex/html/word等。

外部连接性：Maple 的程序可以自动转换为其他语言代码，如

C/Fortran/Java/VB/MATLAB，融合多种开发工具。Maple能够与MATLAB/Simulink, NAG，EXCEl，数据库等工具无缝连接。另外Maple可与CAD系统连接，可通过参数传输完成对CAD模型的数学分析，如统计分析、优化、经验公式计算、公差和单位计算，并自动在CAD系统中完成更新。通过专业工具箱，Maple可与数值计算软件Matlab共享命令、变量等。

二、使用工具的命令，函数

Maple内置5000多个计算命令，深度涵盖广泛的数学和编程主题。

（一）调用函数的命令

调用Maple函数的常用方法一般有两种。分别是maple(statement)和maple('function',arg1,arg2,…)。函数maple用于调用Maple中的符号计算“引擎”和它庞大的函数库，其中的输入参数既可以是符号变量，也可以是数值变量，其输出结构和输入参数的类型相同。

1.maple(statement)

该命令的功能是把对变量、表达式以及函数等的描述（statement）传递到MATLAB的符号运算引擎——Maple V中去，由这个引擎来完成计算任务，返回字符型结果。此函数功能非常强，它可以调用出图像处理函数外的所有函数。例如：

dsolve函数是Maple中解微分方程的函数，可以是单个微分方程，也可以是微分方程组。

【实例1】

>> maple(dsolve('(Dy)^2+y^2=1','y(0)=0.5'))

ans =

[ sin(t+1/6*pi)]

[ cos(t+1/3*pi)]

2. maple('function',arg1,arg2,…)

该命令是MATLAB调用Maple函数库的标准形式。输入参数function 是函数的名，arg1,ar g2,…是function函数的参数，参数个数根据function函数具体确定。

【实例2】

>>maple('mtaylor(exp(x),[x=0],5)')

ans =

mtaylor(exp(x),[x = 0],5)

注意：Maple有些函数在MATLAB启动时，并没有自动调入内存当中。上面的例子中并没有进行分解，所以在调用这些函数时候必须先读取Maple函数库。

【实例3】

>> maple('readlib(mtaylor)')

ans =

proc () local f, k, v, m, n, s, t, w; option `Copyright (c) 1991 by………………

然后我们就重新调用maple('mtaylor(exp(x),[x=0],5)')如下：

>> maple('mtaylor(exp(x),[x=0],5)') ans =

1+x+1/2*x^2+1/6*x^3+1/24*x^4

（二）工具包含的常用函数

三、例题

1.F（x，y); M(x, y) = F(x, y)*(x+y^2); N(x, y) = F(x, y)*xy diff(M(x, y), y) = diff(N(x, y), x)

解：M:=F(x,y)*(x+y^2)

N:=F(x,y)*x*y

eq:=diff(M,y)=diff(N,x)

pdsolve(eq)

F(x,y)=x*G(2x3+3x2y2)，G 是任意函数

2. 求dt

dx =y ，dt dy =-x 的解。 解：S=dsolve('Dx = y ,Dy = -x')

disp([blanks(12),'x',blanks(21),'y']),disp([S.x,S.y])

S =

x: [1x1 sym]

y: [1x1 sym]

x y

[ C1*sin(t)+C2*cos(t), C1*cos(t)-C2*sin(t)]

3. 求解两点边值问题：xy ’’- 3y ’ = 2x ,y(1) = 0 ,y(5) = 0。 解：y =dsolve('x*D2y - 3*Dy = x^2','y(1) = 0,y(5) =0','x')

y =

31/468*x^4-1/3*x^3+125/468

4. 求解微分方程的初值问题dy

dy +3y=8，0x y/ =2的解。 解：p=dsolve('Dy+3*y=8','y(0)=2','x')

p =

8/3-2/3*exp(-3*x)

5.求微分方程 w'(t)=-k*f(t) f'(t)=g*f(t)*[m*w(t)-f(t)]。

建立一个文件test.m

function dx=test(t,x)

dx=zeros(2,1);

dx(1)=k*x(2);

dx(2)=g*x(2)(m*x(1)-x(2));

保存,

在matlab 里输入:

[t,x]=ode45('test',[0 15],[25 2]);

maple 的

restart:with(PDEtools):

alias(w=w(t),f=f(t)):

eq1:=diff(w,t)=k*f;

eq2:=diff(f,t)=g*f*(m*w-f);

dsolve({w,f});

但是方程本身不可积，所以无法求出解。

6. 求解van der Pol 方程：

???-='==+'--''1.0)0(,0)0(

)1(2y y y y y y

> ODE:=diff(y(t),t$2)-(1-y(t)^2)*diff(y(t),t)+y(t)=0;

:= ODE = - + ?? ??????2t ()y t () - 1()y t 2?? ?????()y t ()y t 0

> initvals:=y(0)=0,D(y)(0)=-0.1;

:= initvals , = ()y 00 = ()()D y 0-.1

> F:=dsolve({ODE,initvals},y(t),type=numeric); := F proc () ... end proc

rkf45_x F(0);

??????

??,, = t 0. = ()y t 0. = ??t ()y t -.1 > F(1);

??????

??,, = t 1. = ()y t -.144768589749425608 = ??t ()y t -.178104066128215944

7. 求解：042=+''+'y x y y x

> ODE:=x*diff(y(x),x$2)+diff(y(x),x)+4*x^2*y(x)=0;

:= ODE = + + x ?? ??????2x 2()y x ?? ?????()y x 4x 2()y x 0 > dsolve(ODE,y(x));

= ()y x + _C1?? ????BesselJ ,043x ()/32_C2?? ??

??BesselY ,043x ()/32 > initvals:=y(0)=y0,D(y)(0)=0;

:= initvals , = ()y 0y0 = ()()D y 00

> with(powseries): > sol:=powsolve({ODE,initvals});

:= sol proc () ... end proc powparm

> tpsform(sol,x,16);

- + - + - + y049y0x 3481y0x 6166561y0x 9459049y0x 121613286025

y0x 15()O x 16

8.求解一维谐振子的解：0)(2=-+''y x y ε

> alias(y=y(x)):

> ODE:=diff(y,x$2)+(epsilon-x^2)*y=0;

:= ODE = + ?? ??

????2x 2y () - εx 2y 0 > H:=powsolve(ODE);

:= H proc () ... end proc powparm

> tpsform(H,x,8);

C0C1x 1εC0x 21εC1x 3?? ???? + 1ε2C01C0x 4?? ???? + 1ε2C11C1x 5?? ????- - 1ε?? ???? + 1ε2C01C01εC0 + - - + + + x 6?? ????- - 1ε?? ???? + 1ε2C11C11εC1x 7()O x 8 + +

> H(_k); - - ε()a - _k 2()a - _k 4_k ()

- _k 1 9. 求解物理摆的大幅振动方程：

θθsin g l -=, 其中l 是摆长, θ是摆角, g 是重力加速度。

> ODE:=l*diff(theta(t),t$2)=-g*sin(theta(t));

:= ODE = l ?? ??

????2

t 2()θt -g ()sin ()θt > initvals:=theta(0)=0,D(theta)(0)=v[0]/l; := initvals , = ()θ00 = ()()D θ0v 0

> sol:=dsolve({ODE,initvals},theta(t),type=series);

:= sol = ()θt - + + v 0l t 16g v 0l 2t 31120g v 0() + v 02g l l

4t 5()O t 6 > Order:=11:

> sol:=dsolve({ODE,initvals},theta(t),type=series);

:= sol = ()θt - + - + + v 0

t 1g v 0l 2t 31g v 0() + v 02g l l 4t 51g v 0() + + 11g l v 02g 2l 2v 04l 6t 71g v 0() + + + 57g v 04l 102g 2v 02l 2g 3l 3v 06l 8t 9()O t 11

四、小结

通过对Matlab软件中Maple的学习，我知道了Maple的一些基本功能和用法。