《數(shù)值分析》教案16

《《數(shù)值分析》教案16》由會員上傳分享，免費在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫。

1、3.4常微分方程初值問題數(shù)值解的MATLAB實現(xiàn)在MATLAB中，有多個求解常微分方程數(shù)值解的函數(shù)命令，在此介紹幾個常用的函數(shù)。3.4.1求常微分方程初值問題數(shù)值解的函數(shù)ode23和ode45是求解常微分方程數(shù)值解最常用的兩個函數(shù)，命令中的“ode”是英文常微分方程“OrdinaryDifferentialEquation”的縮寫，它們都采用龍格-庫塔公式進行數(shù)值求解，23和45分別表示使用是2/3階和4/5階龍格-庫塔公式。它們的調(diào)用格式基本相同。在此僅以ode23為例來說明函數(shù)的用法。函數(shù)ode23的調(diào)用格式為：[x,y]=o

2、de23('Fun',Tspan,y0,options)1.該命令適用于一階常微分方程組，，如遇到高階常微分方程，必須先把它們變成一階常微分方程組，即狀態(tài)方程，方可使用。2.輸入?yún)?shù)“Fun”為定義微分方程組的M-函數(shù)文件名，可以在文件名前加寫@，或用英文格式單引號界定文件名。3.在編輯調(diào)試窗口中編寫一階常微分方程組的M-函數(shù)文件時，每個微分方程的格式必須都與一致，即等號左端為待求函數(shù)的一階導(dǎo)數(shù)，右端函數(shù)的變量嚴(yán)格以“先自變量、后函數(shù)”的固定順序輸入，的下標(biāo)表示微分方程的序數(shù)。4.輸入?yún)?shù)“Tspan”規(guī)定了常微分方程的自變量取值

3、范圍，它以矩陣的形式輸入，表示自變量。5.輸入?yún)?shù)“”表示初始條件向量，。微分方程組中的方程個數(shù)必須等于初始條件數(shù)，這是求常微分方程特解所必須的條件。6.輸入?yún)?shù)“option”表示選項參數(shù)，它可由ODESET函數(shù)設(shè)置，較為復(fù)雜。7.輸出參數(shù)為微分方程組解函數(shù)的列表（x和y都是列矩陣），它包含向量x各節(jié)點和與對應(yīng)向量y的第i個分量值(即第i個方程解)，i表示節(jié)點序列數(shù)。表3-1微分方程初值問題的計算函數(shù)函數(shù)名階的精確度如何使用ode45中等首先嘗試此函數(shù)ode23低低精度容差或適度剛性問題ode113從低至高高精度容差ode15s

4、從低至中如果用ode45求解很慢ode23s低低精度容差的線性剛性系統(tǒng)ode23t低問題是適度剛性的ode23tb低用低精度容差求解剛性系統(tǒng)8.輸出參數(shù)缺省時輸出解函數(shù)的曲線，即函數(shù)及其各階導(dǎo)數(shù)的曲線。求解微分方程的命令還有ode45，ode113，ode15s，ode23s，ode23t，ode23tb（見表3-1）等。剛性方程組又稱為Stiff方程組，其特點是出現(xiàn)解的分量級差別很大，給數(shù)值計算帶來很大的困難。在化學(xué)反應(yīng)、電子網(wǎng)絡(luò)和自動控制等領(lǐng)域中經(jīng)常遇到。3.4.2ode23與ode45使用方法舉例【例3-25】求二階常微分方

5、程的通解與特解。解：將二階常微分方程變換成兩個一階常微分方程組：第一種解法：先建立M-文件。%定義輸入,輸出變量和函數(shù)文件名Functiondy=myfun_1(x,y)%明確dy的維數(shù),用微分方程組時不可缺省dy=zeros(2,1);%dy(m)表示y的m階導(dǎo)數(shù);y(n)表示y的第n列dy(1)=y(2);%與方程組中第二個微分方程對應(yīng)dy(2)=2*y(2)-2*y(1)+5*exp(2*x)*sin(x);在命令窗口鍵入：>>[xy]=ode23(@myfun_1,[0,1],[-2;-3])回車得到：x=00.0533…

6、…1.0000y=-2.0000-3.0000-2.1627-3.0957……-1.769112.8927顯示的結(jié)果是自變量x和兩個待求函數(shù)和的對應(yīng)數(shù)據(jù)。其中用“……”代替了許多輸出數(shù)據(jù)。再鍵入：>>x=size(x),y=size(y)x=151y=152表明自變量被分為15個節(jié)點，并計算出了對應(yīng)點上的函數(shù)及其一階導(dǎo)數(shù)的取值。第二種解法：使用內(nèi)聯(lián)函數(shù)inline。在命令窗口鍵入：>>myfun_2=inline('[y(2);2*y(2)-2*y(1)+5*exp(2*x)*sin(x)]','x','y');>>[xy]=od

7、e23(myfun_2,[01],[-2-3])回車得到與第一種方法相同的結(jié)果。在編輯窗口再鍵入（不寫輸出參量）：ode23(@myfun_1,[0,1],[-2;-3])legend('特解函數(shù)','一階導(dǎo)函數(shù)')text(0.8,-2,'特解函數(shù)','FontSize',9)text(0.8,3.5,'一階導(dǎo)函數(shù)','FontSize',9)運行得到圖3-1所示的微分方程的圖示解。圖3-1特解函數(shù)及其導(dǎo)數(shù)曲線寫有輸出參量[xy]時，得出微分方程的數(shù)值解，不寫輸出參量[xy]時，得出的微分方程的圖示解?！纠?-26】解微分方程解

8、：用第一種方法解微分方程，建立M-文件：%定義輸入,輸出變量和函數(shù)名functiondy=myfun_2(x,y)%明確dy的維數(shù),用微分方程組時不可缺省dy=zeros(3,1);%dy(m)表示y的m階導(dǎo)數(shù);dy(1)=y(2);%y(n)表示