數學實驗第八章微分方程

8.1引例解:設所求曲線方程為y=y(x),則又y(x)應滿足條件x=1時y=2，因此可解得C=1兩端對x積分，得例1求xOy平面上過點(1,2)的一條曲線方程，使其在任意點處的斜率為2x。故所求曲線方程為y=x2+1。例2列車上以20米/秒的速度行駛，當制動時，列車獲得加速度-0.4米/秒2，問開始制動后多長時間能停住，以及列車在這段時間行駛了多少路程？等式兩邊同時積分，得解:設列車開始制動后t秒內行駛了s米，則再積分得利用條件t=0,s=0,ds/dt=20，可以解得C1=20,C2=0,因此利用條件列車停住時速度為零，由解得:s=500(米)代入及C1=20可得:t=50(秒)微分方程:例如8.2微分方程的基本概念含有未知函數的導數或微分的方程常微分方程(ODE)(OrdinaryDifferentialEquations)僅含一個自變量的微分方程.微分方程的階:

微分方程中出現的未知函數的最高階導數的階數.線性與非線性微分方程:方程中關于未知函數及其各階導數均是一次的,則稱為線性微分方程.線性微分方程非線性微分方程微分方程的解:代入微分方程能使方程成為恒等式的函數稱為微分方程的解。微分方程解的分類：(1)通解:微分方程的解中含有任意常數,且任意常數的個數與微分方程的階數相同。(2)特解:確定了通解中任意常數以后的解。解的圖象(積分曲線):

微分方程的解曲線。初始條件:

用來確定任意常數的條件。考慮引例1，對于微分方程函數y=x2+C是其通解，y=x2+1是其滿足初始條件x0=1，y(x0)=2的特解。微分方程的(部分)積分曲線(0≤C≤2)見右圖，其中紅色曲線為滿足初始條件的特解曲線。考慮引例2，對于微分方程函數s=-0.2t2+C1t+C2是其通解,s=-0.2t2+20t是其滿足條件t0=0，s(t0)=0,ds/dt(t0)=20的特解.滿足初始條件的積分(特解)曲線見下圖。路程曲線速度曲線

微分方程是研究函數變化規律的有力工具，在科技、工程技術、經濟管理以及生態、環境、人口、交通等各個領域中有著十分廣泛的應用。

很多實際問題的數學模型都可以用微分方程來表示，但建立微分方程只是解決問題的第一步，通常在此基礎上還需要求出微分方程的解并對實際問題進行解釋和加以檢驗。

求解微分方程一般可分為求微分方程的解析解和數值解。微分方程的解

解析解可以用一個確切的代數表達式來表示的解。顯然微分方程的解析解對微分方程的分析和應用都是很方便的。但是存在解析解的微分方程僅僅局限于一些特殊類型，而絕大多數微分方程都是求不出解析解的。

數值解求解不存在解析解的微分方程時必須借助數值解法去求出微分方程的數值解(近似解)。因此，數值解法是求解微分方程的一個十分重要的方法。ODE的解很少能用初等函數及其不定積分的組合表示。例如：方程不能表示為初等函數，故得不到精確解的解是難以求積方程8.3常微分方程數值解法介紹(一)常微分方程數值解的定義在生產和科研中所處理的微分方程往往很復雜且大多得不出一般解析解。而在實際上，對初值問題，一般是要求得到解在若干個點上滿足規定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達式。因此，研究常微分方程的數值解法是十分必要的。對一階微分方程

dy/dx=f(x,y),y(x0)=y0主要方法(1)歐拉方法(向前歐拉公式，向后歐拉公式，梯形公式，改進的歐拉公式)(2)龍格-庫塔方法(二階，三階，四階，五階)在自變量x的取值點列{xn}，根據一定的原理(主要是用差商代替導數)，利用迭代方法求出y(xn)的近似值yn(二)

建立數值解法的一些途徑1、用差商代替導數若步長h較小，則有故有公式：此即歐拉法。設xi+1-xi=h,(i=0,1,…,n-1),可以用以下離散化方法求解微分方程dy/dx=f(x,y),y(x0)=y02、使用數值積分對方程y'=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：故有公式：實際應用時，與歐拉公式結合使用：此即改進的歐拉法。3、使用泰勒公式以此方法為基礎，有龍格-庫塔法、線性多步法等方法。4、數值公式的精度當一個數值公式的截斷誤差y(xi)-yi可表示為O(hk+1)時(k為正整數，h為步長)，稱它是一個k階公式。k越大，則數值公式的精度越高。歐拉法是一階公式，改進的歐拉法是二階公式。龍格-庫塔法有二至五階公式。線性多步法有四階外插公式和內插公式。8.4微分方程求解的MATLAB實現

一、求微分方程(組)的解析解求微分方程(組)解析解的MATLAB命令dsolve

1．求微分方程的通解格式：dsolve('eqn','v')

其中'eqn'

是微分方程表達式構成的字符串，'v'是微分方程表達式中自變量符號。'v'可以省略，省略時默認t為自變量。

在表達微分方程時,用字母D表示微分,D2,D3表示二,三階微分。即在不指定自變量(默認t為自變量)時，Dy=dy/dt，D2y=d2y/dt2,…例8.1求微分方程dy/dx=1+y2的通解

dsolve('Dy=1+y^2')

ans=tan(t-C1)其中t為默認的自變量。如果希望以x為自變量，則需用命令

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

ans=tan(x-C1)注意：該命令中自變量的取法只要不與因變量相同即可。C1cos(x)+C2sin(x)---------------------1/2x例8.2求微分方程

x2y″+xy′+(x2-1/4)y=0的通解s=dsolve('x^2*D2y+x*Dy+(x^2-1/4)*y=0','x')

s=(C1*cos(x)+C2*sin(x))/x^(1/2)注意(1)若方程是右端等于0的形式，可以省略。(2)若不加自變量x，則將把x作為常數求解。(3)所得結果可以用命令pretty進行簡化。例pretty(s)2.求微分方程的特解格式:dsolve('eqn','y(a)=b0,Dy(a)=b1,···','v')其中'eqn'是微分方程表達式構成的字符串，'y(a)=b0,Dy(a)=b1,···'

是微分方程的初始條件。'v'是微分方程表達式中自變量符號。例8.3求微分方程y''+4y'+29y=0滿足條件y(0)=0,y’(0)=15的特解。dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')

ans=3*exp(-2*x)*sin(5*x)3.求解微分方程組

格式：dsolve('eqn1','eqn2',…)例8.4求下面微分方程組的通解。

dx/dt=3x+4y，dy/dt=-4x+3y

[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')x=exp(3*t)*cos(4*t)*C1+exp(3*t)*sin(4*t)*C2y=-exp(3*t)*sin(4*t)*C1+exp(3*t)*cos(4*t)*C2注意：也可以用此命令求滿足條件的特解

[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=1,y(0)=0')x=exp(3*t)*cos(4*t)y=-exp(3*t)*sin(4*t)二、求微分方程(組)的數值解注意求微分方程(組)數值解的MATLAB命令ode微分方程的數值解命令ode僅僅適用于一階微分方程(組)。對高階微分方程必須先將其化為一階微分方程(組)再求解。高階微分方程化為一階微分方程組的方法

對任意高階微分方程

x(n)=f(t,x,x',···,x(n-1))

令y1=xy2=y1'=x'y3=y2'=x''···yn=y'n-1=x(n-1)y1'=y2y2'=y3···yn-1'=ynyn'

=f(t,y1,y2,···,yn)

所以下面只討論一階微分方程組的數值解

一階微分方程組的向量形式表示

令一階微分方程組通常可用向量形式表示為dy/dt=f(t,y)

一階微分方程組的向量形式表示

dy/dt=f(t,y)一階線性微分方程組的向量形式表示

dy/dt=Ay其中求一階微分方程組數值解的MATLAB命令ode

ode45,ode23,ode113,ode15s,ode23s五個命令的格式完全一致,常用的是前兩個。格式：[T,Y]=ode45('fun',Ts,Y0)'fun'：由微分方程組寫成的函數M文件的文件名；Ts=[t1,t2]：方程組解的取值區間；Y0

：方程組解的初始值(n維向量)輸出T是m維列向量(取值從t1到t2)；Y是m×n維矩陣，其第i列是與T相應的yi的值。向量微分方程組dy/dt=f(t,y)的ode輸出[T,Y]=ode45('fun',Ts,Y0,options)

例8.5求解線性微分方程組X'=AX，其中x1'

=x2x2'=4x1+3x2-4x3x3'=x1+2x2+x3

解：建立函數M文件exam1：functiondx=exam1(t,x)dx=zeros(3,1);dx(1)=x(2);dx(2)=4*x(1)+3*x(2)-4*x(3);dx(3)=x(1)+2*x(2)+x(3);取初始條件x(0)=[1;2;3],Ts=[0,20],利用ODE45命令求解如下[T,Y]=ode45('exam1',[0,20],[1;2;3]);利用命令plot(T,Y(:,1))可得積分曲線.例8.6求解Vander

Pol方程x''+(x2-1)x'+x=0

解：令y1=x,y2=y1'=x',則方程化為方程組

y1'

=y2,y2'

=(1-y12)y2-y1建立函數M文件vdp：function

dy=vdp(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=(1-y(1)^2)*y(2)-y(1);取初始條件y(0)=[0;1],Ts=[0,20],利用ODE45命令求解如下Ts=[0,20];y0=[0;1];[T,Y]=ode45('vdp',Ts,y0);積分曲線

由y1=x可知，原方程的解是輸出矩陣Y的第一列Vander

Pol方程x''+(x2-1)x'+x=0滿足