matlab在常微分方程數(shù)值解中的應(yīng)用

1、MATLAB在常微分方程數(shù)值解中的應(yīng)用【摘 要】許多現(xiàn)實問題都可以通過微分方程的形式進(jìn)行表示，傳統(tǒng)解微分方程的方法有近似分析解法、表解法和圖解法，這些方法需對其進(jìn)行大量的假設(shè)，而使得數(shù)學(xué)模型有一定的失真，有一定的局限性。數(shù)值解法利用計算機(jī)，使得求解更精確、效率更高，而MATLAB是一種數(shù)學(xué)軟件包，有高級編程格式，使得計算結(jié)果更具有可信性，因此微分方程的求解及MATLAB在其中的應(yīng)用具有實際意義。本文對常微分方程數(shù)值解問題作進(jìn)一步探討，并應(yīng)用MATLAB對其中難解的改進(jìn)Euler法和Runge-Kutta法進(jìn)行編程實現(xiàn)，程序簡潔、直觀，求解速度快、方法實用性較強(qiáng)?！娟P(guān) 鍵 詞】常微分方程 數(shù)值解

2、 MATLAB Euler法 龍格-庫塔方法 ode45 ode15sMatlab in ordinary differential equation numerical solution of applicationYang Hua Zhang Lei【Abstract】Many practical problems can be using differential equations in the form of representation, the traditional method of solving differential equations are similar ana

3、lysis method, table method and graphical method, these methods to carry on the large amounts of hypothesis, so that the mathematical model has certain distortion, have certain limitation. Numerical solution of using a computer, make solving more accurate and more efficient, and MATLAB is a kind of m

4、athematics software package, with advanced programming format, making calculation result is more credibility, therefore differential equation and solution of the MATLAB in one of the application of practical significance. This paper numerical solution of differential equation problem further discuss

5、ion, and the application of MATLAB in which the difficult solution improvement Euler method and Runge - Kutta method on the programming, the program is concise, intuitive and solution speed, method of practical stronger.【Key words】ordinary differential equation，numerical solution，Matlab，Euler method

6、，Runge-Kutta method 【引 言】微分方程的概念：未知的函數(shù)以及它的某些階的導(dǎo)數(shù)連同自變量都由一已知方程聯(lián)系在一起的方程稱為微分方程。如果未知函數(shù)是一元函數(shù)，稱為常微分方程。常微分方程的一般形式為微分方程是數(shù)學(xué)科學(xué)聯(lián)系實際問題的主要橋梁之一，它是含有未知數(shù)及其導(dǎo)數(shù)的方程。常微分方程的求解是現(xiàn)代科學(xué)研究和工程技術(shù)中經(jīng)常遇到的問題，然而，從實際問題中建立起來的微分方程往往具有非常復(fù)雜的形式，有寫解析式難以計算，有的則根本不能用解析式來表達(dá)。在實際上對初值問題，一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值，或者得到一個滿足精確度要求的便于計算的表達(dá)式，所以利用數(shù)值解法求解實際問

7、題顯得非常重要。一階常微分方程的初值，其一般形式為dydx=f（x，y）ya=y0（axb） （1）在本文的討論中，總假設(shè)函數(shù)f（x，y）連續(xù)，且關(guān)于y滿足萊布尼茲條件，即存在常數(shù)L，使得fx，y-f（x，y）Ly-y由微分方程理論可知，初值問題式（1）的解必定存在唯一。若給出節(jié)點xn=a+nh （n=0,1,2 ），其中h為步長，設(shè)y（xn）代表方程（1）的精確解在xn的值，yn代表某種算法（忽略計算的舍入誤差）得到的y（xn）近似值，所謂數(shù)值解法，就是尋求y（x）在系列離散點a=b處的近似值yn（n=1,2,），求解過程是順著節(jié)點的次序一步步的向前推進(jìn)，即按遞推公式由已知的y1，y2，y3

8、 yi求出yi+1。1 建立數(shù)值解法的一些途徑1.1、用差商代替導(dǎo)數(shù) 若步長較小，則有固有公式：此即Euler法。1.2、使用數(shù)值積分對微分方程兩邊積分得故有公式： 和yn+1=yn+h(yn-2xn/yn)yn+1=yn+hyn-2xnyn+yn+1-2xnyn+12此即改進(jìn)的Euler法。1.3四階Runge-Kutta法以Euler法為基礎(chǔ)，繼續(xù)推廣和處理，可得一、二、三、四階Runge-Kutta格式，最常用的一種經(jīng)典Runge-Kutta格式的具體形式如下： yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h2,yn+k12)k3=hf(x

9、n+h2,yn+k22)k4=hf(xn+h,yn+k3)(n=0,1,2,)2利用matlab求解一階常微分方程的初值解在MATLAB中，由函數(shù)dsolve()解決常微分方程（組）的求解問題，其具體格式如下：X=dsolve(f1,f2,)函數(shù)dsolve用來解符號常微分方程、方程組，如果沒有初始條件，則求出通解，如果有初始條件，則求出特解。有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，我們需要尋求方程的數(shù)值解，在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，我們將其統(tǒng)稱為solver，其一般格式為：T,Y=solver(odefx,t,y0)該函數(shù)表

10、示在區(qū)間t=t0,tf上，用初始條件y0求解顯式常微分方程。solver為命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，這些命令各有特點。列表說明如下：求解器ODE類型特點說明ode45非剛性一步算法，4,5階Runge-Kutta方法累積截斷誤差大部分場合的首選算法ode23非剛性一步算法，2,3階Runge-Kutta方法累積截斷誤差使用于精度較低的情形ode113非剛性多步法，Adams算法，高低精度均可達(dá)到計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法，Gears反向數(shù)值積分，精度中等若

11、ode45失效時，可嘗試使用ode23s剛性一步法，2階Rosebrock算法，低精度。當(dāng)精度較低時，計算時間比ode15s短odefx為顯式常微分方程中的，t為求解區(qū)間，要獲得問題在其他指定點上的解，則令t=t0,t1,t2,（要求單調(diào)），y0初始條件。3 根據(jù)數(shù)值解法的思想編程實現(xiàn)MATLAB 中有幾個專門用于求解常微分方程的函數(shù),它們的設(shè)計思想基于Runge - Kutta方法,基本設(shè)計思想為：從改進(jìn)的歐拉方法比歐拉方法精度高的緣由著手,如果在區(qū)間 x1 , xi+1 多取幾個點的斜率值,然后求取平均值,則可以構(gòu)造出精度更高的計算方法。 這些函數(shù)主要包括：ode45、ode23、ode1

12、5s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函數(shù)ode45,該函數(shù)采用變步長四階五階Runge - Kutta法求數(shù)值解,并采用自適應(yīng)變步長的求解方法。ode23采用二階三階Runge - Kutta法求數(shù)值解,與ode45類似,只是精度低一些。ode15s用來求剛性方程組。 ode45函數(shù)的調(diào)用格式為: tout, yout = ode45 ( yprime, t0, tf , y0) 其中yprime是表示f ( t, y)的M文件名, t0表示自變量的初始值, tf表示自變量的終值, y0表示初始向量值. 輸出向量tout表示節(jié)點( t0 ; t1 ;

13、 ; tn ) ,輸出矩陣yout表示數(shù)值。例1Rossler方程的數(shù)值解法 Rossler方程: 求解過程如下: (1)建立自定義函數(shù) 用edit命令建立自定義函數(shù)名為rossler. m, 內(nèi)容為: function dx = rossler ( t, x) dx=-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3);(2) 調(diào)用對微分方程數(shù)值解ode45函數(shù)求解 用edit命令建立一個命令文件rossler1. m, 內(nèi)容為： >>x0=0;0;0; >>t,y=ode45(rossler,0,100,x0); plot(t,y

14、); >>figure; >>plot3(y(:,1),y(:,2),y(:,3) 得 如果微分方程由一個或多個高階常微分方程給出，要得到該方程的數(shù)值解,可以將方程轉(zhuǎn)換成一階常微分方程組。假設(shè)高階常微分方程的一般形式為y( n) = f ( t, y, y, , y( n - 1) )，而且函數(shù)y(t)的各階導(dǎo)數(shù)初值為y(0)，y(0) , y( n - 1) (0) 可以選擇一組變量令: x1 = y, x2 = y, xn = y( n - 1) ,我們就可以把原高階常微分方程轉(zhuǎn)換成下面的一階常微分方程組形式： 而且初值x1(0)=y(0),x2(0)=y(0),，

15、xn(0)=(0)。 轉(zhuǎn)換以后就可以求原高階常微分方程的數(shù)值解了。例2 求微分方程，,的數(shù)值解。對方程進(jìn)行變換,選擇變量 (1) 建立自定義函數(shù) 用edit命令建立自定義函數(shù)名為f.m,內(nèi)容為: function y =f(t,x) y=x(2);x(3);-t2*x(2)*x(1)2-t*x(1)*x(3)+exp(-t*x(1);(2)調(diào)用對微分方程數(shù)值解ode45函數(shù)求解 用edit命令建立一個命令文件f2. m,內(nèi)容為: >>x0=2;0;0; >>t,y =ode45(f,0,10,x0);plot(t,y); >>figure; >>plot3(y(:,1),y(:,2), y(:,3) 得 求剛性微分方程,的數(shù)值解。 調(diào)用對微分方程數(shù)值解ode15s函數(shù)求解 >>x0=2;0;0; >>t,y=ode15s(f,0,10,x0);plot(t,y(:,1) >>figure; >>plot(t,y(:,2) 得 x1 ( t) x2 ( t) 用剛性方程求解函數(shù)可以快速求出該方程的數(shù)值解,并且畫出兩個狀態(tài)變量的時間曲線。x1 ( t) 曲線變化比較平滑, x2 ( t) 曲線變化在某些點上較快.