來源：網絡論壇轉載(VB資料庫)

常微分方程求解

微分方程求解是數學研究與應用的一個重點和難點. Maple能夠顯式或隱式地解析地求解許多微分方程求解.

在常微分方程求解器dsolve中使用了一些傳統的技術例如laplace變換和積分因子法等,

函數pdesolve則使用諸如特征根法等經典方法求解偏微分方程. 此外, Maple還提供了可作攝動解的所有工具,

例如Poincare-Lindstedt法和高階多重尺度法.

幫助處理常微分方程(組)的各類函數存于Detools軟件包中, 函數種類主要有：可視化類的函數, 處理寵加萊動態系統的函數,

調整微分方程的函數, 處理積分因子、李對稱法和常微分方程分類的函數, 微分算子的函數, 利用可積性與微分消去的方法簡化微分方程的函數,

以及構造封閉解的函數等.

更重要的是其提供的強大的圖形繪制命令Deplot能夠幫助我們解決一些較為復雜的問題.

2.1 常微分方程的解析解

求解常微分方程最簡單的方法是利用求解函數dsolve. 命令格式為：

dsolve(ODE);

dsolve(ODE, y(x),

extra_args);

dsolve({ODE, ICs}, y(x),

extra_args);

dsolve({sysODE, ICs}, {funcs}, extra_args);

其中, ODE—常微分方程, y(x)—單變量的任意變量函數, Ics—初始條件, {sysODE}—ODE方程組的集合,

{funcs}—變量函數的集合,

extra_args—依賴于要求解的問題類型.

例如, 對于一階常微分方程 可用dsolve直接求得解析解：

> ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x))-y(x);

> dsolve(ODE,y(x));

可以看出, dsolve的第一個參數是待求的微分方程, 第二個參數是未知函數. 需要注意的是, 無論在方程中還是作為第二個參數,

未知函數必須用函數的形式給出(即：必須加括號, 并在其中明確自變量), 這一規定是必須的,

否則Maple將無法區分方程中的函數、自變量和參變量, 這一點和我們平時的書寫習慣不一致. 為了使其與我們的習慣一致,

可用alias將函數用別稱表示：

> alias(y=y(x));

> ODE:=x*diff(y,x)=y*ln(x*y)-y;

> dsolve(ODE,y);

函數dsolve給出的是微分方程的通解,

其中的任意常數是用下劃線起始的內部變量表示的.

在Maple中, 微分方程的解是很容易驗證的,

只需要將解代入到原方程并化簡就可以了.

> subs(%,ODE);

> assume(x,real): assume(_C1,real):

> simplify(%);

> evalb(%);

evalb函數的目的是對一個包含關系型運算符的表達式使用三值邏輯系統求值, 返回的值是true, false和FAIL.

如果無法求值, 則返回一個未求值的表達式. 通常包含關系型運算符“=, <>, ,

>=”的表達式在Maple中看作是代數方程或者不等式. 然而,

作為參數傳遞給evalb或者出現在if或while語句的邏輯表達式中時, 它們會被求值為true或false. 值得注意的是,

evalb不化簡表達式, 因此在使用evalb之前應將表達式化簡, 否則可能會出錯.

再看下面常微分方程的求解：

> alias(y=y(x)):

> ODE:=diff(y,x)=sqrt(y^2+1);

> dsolve(ODE,y);

函數dsolve對于求解含有未知參變量的常微分方程也完全可以勝任:

> alias(y=y(x)):

> ODE:=diff(y,x)=-y/sqrt(a^2-y^2);

> sol:=dsolve(ODE,y);

由此可見, 對于不能表示成顯式結果的微分方程解, Maple盡可能將結果表示成隱式解. 另外, 對于平凡解y=0常常忽略,

這一點應該引起注意.

dsolve對于求解微分方程初值問題也十分方便的:

> ODE:=diff(u(t),t$2)+omega^2*u(t)=0;

> dsolve({ODE,u(0)=u0,D(u)(0)=v0},u(t));

2.2 利用積分變換求解微分方程

對于特殊的微分方程, 我們還可以指定dsolve利用積分變換方法求解,

只需要在dsolve中加入可選參數method=transform即可. 其中transform是積分變換,

可以是laplace、fourier、fouriercos或者fouriersin變換.

作為例子, 我們來看一個具有阻尼的振子在階躍沖擊(Heaviside函數)下的響應：

>

ODE:=diff(u(t),t$2)+2*d*omega*diff(u(t),t)+omega^2*u(t)=Heaviside(t);

> initvals:=(u(0)=u[0],D(u)(0)=v[0]);

>

solution:=dsolve({ODE,initvals},u(t),method=laplace);

Maple給出了問題的通解, 但沒有區分自由振動(d=0)、欠阻尼(01)的情況.

下面加以區分求解：

> assume(omega>0):

> simplify(subs(d=0,solution));

> K:=subs(d=1/5,u[0]=1,v[0]=1,solution);

> with(plots):

>

plot3d(rhs(%%),omega=2/3..4/3,t=0..20,style=hidden,orientation=[-30,45],axes=framed);

對于d=1的情況, 可可用下式獲得結果：

> limit(rhs(solution),d=1);

再如下例：

>

diff(u(t),t$2)+3*diff(u(t),t)+2*u(t)=exp(-abs(t));

> dsolve(%,u(t),method=fourier);

2.3 常微分方程組的求解

函數dsolve不僅可以用來求解單個常微分方程, 也可以求解聯立的常微分方程組. 特別是對于線性微分方程組,

由于數學上具有成熟的理論, Maple的求解也是得心應手. 其命令格式為：

dsolve( {eqn1, eqn2, …, ini_conds}, {vars});

其中, ini_conds是初始條件.

>

eqn1:={diff(x(t),t)=x(t)+y(t),diff(y(t),t)=y(t)-x(t)};

> dsolve(eqn1,{x(t),y(t)});

> eqn2:=2*diff(x(t),t$2)+2*x(t)+y(t)=2*t;

> eqn3:=diff(y(t),t$2)+2*x(t)+y(t)=t^2+1;

> dsolve({eqn2, qn3, x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=0},

{x(t),y(t)} );

2.4 常微分方程的級數解法

1) 泰勒級數解法

當一個常微分方程的解析解難以求得時, 可以用Maple求得方程解的級數近似, 這在大多數情況下是一種非常好的方法.

級數解法是一種半解析半數值的方法.

泰勒級數法的使用命令為：

dsolve({ODE,Ics}, y(x), 'series'); 或dsolve({ODE,Ics}, y(x),

'type=series');

下面求解物理擺的大幅振動方程： , 其中l是擺長,?是擺角,

g是重力加速度.

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

> initvals:=theta(0)=0,D(theta)(0)=v[0]/l;

>

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

> Order:=11:

>

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

2) 冪級數解法

對于一個符號代數系統來說, 冪級數是必不可少的微分方程求解工具. 冪級數求解函數powsolve存于工具包powseries中.

但是, 這一求解函數的使用范圍很有限,

它只可以用來求解多項式系數的線性常微分方程或方程組,其求解命令為：powseries[function]

(prep)或直接載入軟件包后用function(prep),

prep為求解的線性微分方程及其初值.

例：求解：

>

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

> dsolve(ODE,y(x));

> initvals:=y(0)=y0,D(y)(0)=0;

> with(powseries):

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

> tpsform(sol,x,16);

也可以用powsolve給出的函數直接獲得用遞歸形式定義的冪級數系數, 不過參數必須用_k,

這是powsolve使用的臨時變量.

> sol(_k);

例：求解一維諧振子的解：

> alias(y=y(x)):

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

> H:=powsolve(ODE);

> tpsform(H,x,8);

> H(_k);

2.5 常微分方程的數值解法

在對微分方程的解析解失效后, 可以求助于數值方法求解微分方程. 數值求解的好處是只要微分方程的條件足夠多時一般都可求得結果,

然而所得結果是否正確則必須依賴相關數學基礎加以判斷.

調用函數dsolve求常微分方程初值問題的數值解時需加入參數type=numeric.

另一方面, 常微分方程初值問題數值求解還可以選擇算法, 加入參數“method=方法參數”即可,

方法參數主要有：

rkf45：4~5階變步長Runge-Kutta-Fehlberg法

dverk78：7~8階變步長Runge-Kutta-Fehlberg法

classical：經典方法, 包括向前歐拉法, 改進歐拉法, 2、3、4階龍格庫塔法,

Sdams-Bashford方法等

gear：吉爾單步法

mgear：吉爾多步法

2.5.1變步長龍格庫塔法

下面用4~5階Runge-Kutta-Fehlberg法求解van der Pol方程：

>

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

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

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

此時, 函數返回的是一個函數,

可以在給定的數值點上對它求值：

> F(0);

> F(1);

可以看到,

F給出的是一個包括t、y(t)、D(y)(t)在內的有序表, 它對于每一個時間點可以給出一組數值表達式. 有序表的每一項是一個等式,

可對其作圖描述.

> plot('rhs(F(t)[2])', t=0..15, title="solution of the Van de

Pol's Equation");

> plots[odeplot](F,[t,y(t)],0..15,title="solution of the Van de

Pol's Equation");

2.5.2吉爾法求解剛性方程

在科學和工程計算中, 常常會遇到這樣一類常微分方程問題, 它可以表示成方程組： , 稱其為剛性方程, 其解的分量數量相差很大,

分量的變化速度也相差很大. 如果用常規方法求解, 為了使變量有足夠高的精度, 必須取很小的步長, 而為了使慢變分量達到近似的穩態解,

則需要很長的時間, 這樣用小步長大時間跨度的計算, 必定造成龐大的計算量, 而且會使誤差不斷積累.

吉爾法是專門用來求解剛性方程的一種數值方法.

>

ODE:=diff(u(t),t)=-2000*u(t)+999.75*v(t)+1000.25,diff(v(t),t)=u(t)-v(t);

> initvals:=u(0)=0,v(0)=-2;

>

ansl:=dsolve({ODE,initvals},{u(t),v(t)},type=numeric,method=gear);

> ansl(10,0);

> p1:=plots[odeplot]

(ansl,[t,u(t)],0..20,color=red):

p2:=plots[odeplot]

(ansl,[t,v(t)],0..20,color=blue):

plots[display] ({p1,p2}, title="Solution of a stiff

equation");

2.5.3 經典數值方法

Maple中常微分方程數值解法中有一類被稱作是“經典”(classical)方法. 當然, 稱其為經典方法不是因為它們常用或是精度高,

而是因為它們的形式簡單, 經常被用于計算方法課上的教學內容. 它們是一些常見的固定步長方法,

在dsolve中用參數method=classical[方法名稱], 如果不特別指出, 將默認采用向前歐拉法.

主要有：

foreuler：向前歐拉法(默認)

hunform：Heun公式法(梯形方法, 改進歐拉法)

imply：改進多項式法

rk2：二階龍格庫塔法

rk3：三階龍格庫塔法

rk4：四階龍格庫塔法

adambash：Adams-Bashford方法(預測法)

abmoulton：Adams-Bashford-Moulton方法(預測法)

下面給出微分方程數值方法的參數表：

參數名

參數類型

參數用途

參數用法

initial

浮點數的一維數組

指定初值向量

number

正整數

指定向量個數

output

'procedurelist'(默認)

或'listprocedure'

指定生成單個函數或多個函數的有序表

Procedurelis：單個函數,返回有序表

Listprocedure：函數的有序表

procedure

子程序名

用子程序形式指定第一尖常微分方程組的右邊部分

參數1：未知函數的個數

參數2：自變量

參數3：函數向量

參數4：導函數向量

start

浮點數

自變量起始值

startinit

布爾量(默認FALSE)

指定數值積分是否總是從起始值開始

對dverk78不適用

value

浮點數向量(一維數組)

指定需要輸出函數值的自變量數值點

如果給定,結果是一個的矩陣.元素[1,1]是一個向量,含自變量名和函數名稱;元素[2,1]是一個數值矩陣,其中第一列value的輸入相同,其他列中是相應的函數值

另外, 還有一些特殊的附加參數：

maxfun：整數類型, 用于最大的函數值數量, 默認值50000, 為負數時表示無限制

corrections：正整數類型, 指定每步修正值數量, 在abmoulton中使用, 建議值≤4

stepsize：浮點數值, 指定步長

下面看一個簡單的例子：

> ODE:=diff(y(x),x)=y(x)-2*x/y(x);

> initvals:=y(0)=1;

>

sol1:=dsolve({ODE,initvals},y(x),numeric,method=classical,stepsize=0.1,start=0);

而其解析解為：

> sol2:=dsolve({diff(y(x),x)=y(x)-2*x/y(x), y(0)=1},

y(x));

將兩者圖形同時繪制在同一坐標系中比較, 可以發現,

在經過一段時間后, 歐拉法的數值結果會產生較大誤差.

> plot({rhs(sol2),'rhs(sol1(x)[2])'},x=0..2);

求解微方程, 無論使用什么方法或者加入什么選項,

求解完成后必須利用相關數學知識進行邏輯判斷, 絕對不對簡單迷信Maple給出的結果, 否則很有可能得到一個對于方程本身也許還看得過去,

但在數學或者物理意義上不合理的解.

2.6攝動法求解常微分方程

由于微分方程求解的復雜性, 一般微分方程常常不能求得精確解析解, 需要借助其它方法求得近似解或數值解, 或者兩種方法兼而有之.

攝動法是重要的近似求解方法.

攝動法又稱小參數法, 它處理含小參數 的系統, 一般當 =0時可求得解x0. 于是可把原系統的解展成 的冪級數 ,

若這個級數當?0時一致收斂，則稱正則攝動, 否則稱奇異攝動.

攝動法的種類繁多, 最有代表性的是龐加萊—林斯泰特(Poicare-Lindstedt)法, 在此, 我們以該方法求解van der

Pol方程：

當 =0時該方程退化為數學單擺的常微分方程, 當 =1時為3.5討論的情況, 對任意 , 該微分方程擁有一個漸進穩定的周期解,

稱為極限環.

由于van der Pol方程中沒有顯式的時間依賴項, 不失一般性, 設初值為y(0)=0. 在龐加萊—林斯泰特法中,

時間通過變換拉伸：

, 其中

對于 , van der Pol方程變為：

restart:

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

>

ODE:=DEtools[Dchangevar]({t=tau/omega,y(t)=y(tau)},%,t,tau);

> e_order:=6:

> macro(e=epsilon,t=tau):

> alias(seq(y=eta(tau),i=0..e_order)):

> e:=()->e:

> for i from 0 to e_order do

eta:=t->eta(t)

od:

> omega:=1+sum('w*e^i','i'=1..e_order);

> y:=sum('eta*e^i','i'=0..e_order);

> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}):

> for i from 0 to e_order do

ode:=coeff(lhs(deqn),e,i)=0

od:

> ode[0];

> ode[1];

> ode[2];

>

dsolve({ode[0],eta[0](0)=0,D(eta[0])(0)=C[1]},eta[0](t));

> eta[0]:=unapply(rhs(%),t);

> ode[1];

> map(combine,ode[1],'trig');

> ode[1]:=map(collect,%,[sin(t),cos(t)]);

>

dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);

> map(collect,%,[sin(t),cos(t),t]);

>

solve({coeff(lhs(ode[1]),sin(t))=0,coeff(lhs(ode[1]),cos(t))=0});

> w[1]:=0:C[1]:=-2:

> ode[1];

>

dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);

> eta[1]:=unapply(rhs(%),tau);

> map(combine,ode[2],'trig'):

>

ode[2]:=map(collect,%,[sin(t),sin(3*t),cos(t),cos(3*t)]);

>

solve({coeff(lhs(ode[2]),sin(t))=0,coeff(lhs(ode[2]),cos(t))=0});

> assign(%):

>

dsolve({ode[2],eta[2](0)=0,D(eta[2])(0)=C[3]},eta[2](t),method=laplace):

> eta[2]:=unapply(rhs(%),t);

> for i from 0 to e_order do

map(combine,ode,'trig'):

ode:=map(collect,%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):

solve({coeff(lhs(ode),sin(t))=0,coeff(lhs(ode),cos(t))=0}):

assign(%):

dsolve({ode,eta(0)=0,D(eta)(0)=C[i+1]},eta(t),method=laplace):

collect(%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):

eta:=unapply(rhs(%),t);

od:

> omega;

> y(t):

> y:=unapply(simplify(y(t),{e^e_order=0}),t):

> e:=1:y(t);

> plot(y(t),t=0..15);

在該問題的求解過程中, 前半部分我們按照交互式命令方式輸入, 也就是把數學邏輯推理的過程“翻譯”成Maple函數, 而在后半部分,

則采用程序設計方式書寫了數學推導過程, 這是應用Maple解決實際問題的兩種方式. 前一種方法只需了解Maple函數即可應用,

而后一種程序設計方式則需掌握Maple程序設計語言. 但是, 不論是那一種方式, 數學基礎總是最重要的.

總結

以上是生活随笔為你收集整理的maple 解代数方程组得多项式_Maple笔记2--常微分方程求解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。