版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請(qǐng)進(jìn)行舉報(bào)或認(rèn)領(lǐng)
文檔簡介
1、求微分方程的解數(shù)學(xué)實(shí)驗(yàn)1;.2q 自牛頓發(fā)明微積分以來,微分方程在描述事物運(yùn)動(dòng)規(guī)律上已發(fā)揮了重要的作用。自牛頓發(fā)明微積分以來,微分方程在描述事物運(yùn)動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問題通過數(shù)學(xué)建模所得到的方程,絕大多數(shù)是微分方程。實(shí)際應(yīng)用問題通過數(shù)學(xué)建模所得到的方程,絕大多數(shù)是微分方程。q 由于實(shí)際應(yīng)用的需要,人們必須求解微分方程。然而能夠求得解析解的微分方程由于實(shí)際應(yīng)用的需要,人們必須求解微分方程。然而能夠求得解析解的微分方程十分有限,絕大多數(shù)微分方程需要利用數(shù)值方法來近似求解。十分有限,絕大多數(shù)微分方程需要利用數(shù)值方法來近似求解。q 本實(shí)驗(yàn)主要研究如何用本實(shí)驗(yàn)主要研究如何用 Matlab
2、 來計(jì)算微分方程(組)的數(shù)值解,并重點(diǎn)介紹一個(gè)來計(jì)算微分方程(組)的數(shù)值解,并重點(diǎn)介紹一個(gè)求解微分方程的基本數(shù)值解法求解微分方程的基本數(shù)值解法Euler折線法。折線法。問題背景和實(shí)驗(yàn)?zāi)康膯栴}背景和實(shí)驗(yàn)?zāi)康?q 考慮一維經(jīng)典初值問題考慮一維經(jīng)典初值問題00 , , ( ,)() , dyf x yy xyxa bdx u 基本思想:用差商代替微商基本思想:用差商代替微商根據(jù)根據(jù) Talyor 公式,公式,y(x) 在點(diǎn)在點(diǎn) xk 處有處有2()()() ()()kkky xy xxxyxOx 1kkhxx 11()()()()( )kkkkky xy xy xy xdyO hdxhhx 21()
3、()()()kkky xy xhyxO h Euler 折線法折線法4初值問題的初值問題的Euler折線法折線法q 具體步驟:具體步驟:等距剖分:等距剖分:0121nnaxxxxxb 步長:步長:1 0 1 21()/, , ,kkhxxknnba u 分割求解區(qū)間分割求解區(qū)間u 差商代替微商差商代替微商1()()kkky xy xdydxhx 1 ()()()kkky xy xh yx 得方程組:得方程組:0011 ()(,)kkkkkkyy xyyh f xyxxh 分割求解區(qū)間,差商代替微商,解代數(shù)方程分割求解區(qū)間,差商代替微商,解代數(shù)方程 為分割點(diǎn)為分割點(diǎn) 0nkkx k = 0, 1
4、, 2, ., n-1yk 是是 y (xk) 的近似的近似5Euler 折線法舉例折線法舉例例:用例:用 Euler 法解初值問題法解初值問題22 0 201 , ( )dyxydxyxy 取步長取步長 h = (2 - 0)/n = 2/n,得差分方程得差分方程00110 1 2 ,(,)(/)kkkkkkkkkkxyyyh f xyyh yxyxxh 當(dāng)當(dāng) h=0.4,即即 n=5 時(shí),時(shí),Matlab 源程序見源程序見 fulu1.m解:解:6Euler 折線法源程序折線法源程序clearf=sym(y+2*x/y2);a=0; b=2;h=0.4;n=(b-a)/h+1; % n=(
5、b-a)/h;x=0; y=1;szj=x,y;for i=1:n-1 % i=1:n y=y+h*subs(f,x,y,x,y);%第一次代初值第一次代初值 x=x+h; szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),or-)7 Euler折線法舉例(續(xù))折線法舉例(續(xù))解析解:解析解:1 3352233/xyex 解析解解析解近似解近似解8Runge-Kutta 方法方法q 為了減小誤差,可采用以下方法:為了減小誤差,可采用以下方法:u 讓步長讓步長 h 取得更小一些;取得更小一些;u 改用具有較高精度的數(shù)值方法:改用具有較高精度的數(shù)值方法:q 龍格龍格
6、-庫塔方法庫塔方法Runge-Kutta (龍格龍格-庫塔庫塔) 方法方法u 是一類求解常微分方程的數(shù)值方法是一類求解常微分方程的數(shù)值方法u 有多種不同的迭代格式有多種不同的迭代格式9Runge-Kutta 方法方法q 用得較多的是用得較多的是 四階四階R-K方法(教材第方法(教材第 79 頁)頁)00111234 (22)/6(),kkkkyy xxxhyyh LLLL 12132432222(,)(/ ,/ )(/ ,/ )(,)kkkkkkkkLf xyLf xhyhLLf xhyhLLf xh yhL 其中其中10四階四階 R-K 方法源程序方法源程序clear;f=sym(y+2*x
7、/y2);a=0; b=2; h=0.4;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1; szj=x,y;for i=1:n-1 % i=1:n l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x,y;endplot(szj(:,1),szj(:,2), dg-)11Runge-Kutta 方法方法12Euler 法與法與
8、 R-K法誤差比較法誤差比較13Matlab 解初值問題解初值問題q 用用 Maltab自帶函數(shù)自帶函數(shù) 解初值問題解初值問題u 求解析解:求解析解:dsolveu 求數(shù)值解:求數(shù)值解: ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb14dsolve 求求解析解解析解q dsolve 的使用的使用y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中其中 y 為輸出,為輸出, eq1、eq2、.為微分方程,為微分方程,cond1、cond2、.為初值條件,為初值條件,v 為自變量。為自變量。例例 1:求微分方程:求
9、微分方程 的通解,并驗(yàn)證。的通解,并驗(yàn)證。22xdyxyxedx y=dsolve(Dy+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2)y = 1/2*exp(-x2)*x2+exp(-x2)*C115dsolve 的使用的使用q 幾點(diǎn)說明幾點(diǎn)說明l 如果省略初值條件,則表示求通解;如果省略初值條件,則表示求通解;l 如果省略自變量,則默認(rèn)自變量為如果省略自變量,則默認(rèn)自變量為 t dsolve(Dy=2*x,x); dy/dx = 2xdsolve(Dy=2*x); dy/dt = 2xl 若找不到解析解,則返回其積分形式。若找不
10、到解析解,則返回其積分形式。l 微分方程中用微分方程中用 D 表示對(duì)表示對(duì) 自變量自變量 的導(dǎo)數(shù),如:的導(dǎo)數(shù),如:Dy y; D2y y; D3y y16dsolve 舉例舉例例例 2:求微分方程求微分方程 在初值條件在初值條件 下的特解,并畫出解下的特解,并畫出解函數(shù)的圖形。函數(shù)的圖形。0 xxyye y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);12( )ye 17dsolve 舉例舉例例例3:求微分方程組:求微分方程組 在初值條件在初值條件 下的特解,并畫出解函數(shù)的圖形。下的特解,并畫出解函數(shù)的圖形。530tdxxyedtdyxy
11、dt x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t)ezplot(x,y,0,1.3);0010|ttxy 18v dsolve 舉例舉例例:例:x,y=dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t) r = dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t)這里返回的這里返回的 r 是一個(gè)是一個(gè) 結(jié)構(gòu)類型結(jié)構(gòu)類型 的數(shù)據(jù)的數(shù)據(jù)r.x %查看解函數(shù)查看解函數(shù) x(t)r.y %查看解函數(shù)查看解函數(shù) y(t)只有很少一部分微分方程(組)能求出
12、解析解。只有很少一部分微分方程(組)能求出解析解。大部分微分方程(組)只能利用大部分微分方程(組)只能利用數(shù)值方法數(shù)值方法求數(shù)值解。求數(shù)值解。dsolve的輸出個(gè)數(shù)只能為一個(gè)的輸出個(gè)數(shù)只能為一個(gè) 或或 與方程個(gè)數(shù)相等。與方程個(gè)數(shù)相等。19Matlab函數(shù)函數(shù)數(shù)值求解數(shù)值求解T,Y = solver(odefun,tspan,y0)其中其中 y0 為初值條件,為初值條件,tspan為求解區(qū)間;為求解區(qū)間;Matlab在數(shù)值求解時(shí)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割,在數(shù)值求解時(shí)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割,T (向量向量) 中返回的是分割點(diǎn)的值中返回的是分割點(diǎn)的值(自變量自變量),Y (向量向量) 中返回的是解函數(shù)
13、在這些分割點(diǎn)上的函數(shù)中返回的是解函數(shù)在這些分割點(diǎn)上的函數(shù)值。值。solver 為為Matlab的的ODE求解器(可以是求解器(可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb)沒有一種算法可以有效地解決所有的沒有一種算法可以有效地解決所有的 ODE 問題,因此問題,因此MATLAB 提供了多種提供了多種ODE求求解器解器,對(duì)于不同的對(duì)于不同的ODE,可以調(diào)用不同的,可以調(diào)用不同的求解器求解器。20參數(shù)說明參數(shù)說明odefun 為顯式常微分方程,可以用命令為顯式常微分方程,可以用命令 inline 定義,或在函數(shù)文件中定義,然后通過定義,或在
14、函數(shù)文件中定義,然后通過函數(shù)句柄調(diào)用。函數(shù)句柄調(diào)用。fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);注:也可以在注:也可以在 tspan 中指定對(duì)求解區(qū)間的分割,如:中指定對(duì)求解區(qū)間的分割,如:x,y=ode23(fun,0:0.1:0.5,1); %此時(shí)此時(shí) x=0:0.1:0.5T,Y = solver(odefun,tspan,y0) 求初值問題求初值問題 的數(shù)值解,求解范圍為的數(shù)值解,求解范圍為 0,0.5222201( )dyyxxdxy 例例 4:21數(shù)值求解舉例數(shù)值求解舉例如果需求解的問題是高階常微分方程,則需將其化為一階
15、常微分方程組,此時(shí)需用函數(shù)文如果需求解的問題是高階常微分方程,則需將其化為一階常微分方程組,此時(shí)需用函數(shù)文件來定義該常微分方程組。件來定義該常微分方程組。122212112101 00 7/()( ),( ),dxdtxdxdtxxxxx 令令 ,則原方程可化為,則原方程可化為12,dyxy xdt 求解求解 Ver der Pol 初值問題初值問題2221001 00 7()( ),( ),d ydyyydtdtyy 例例 5:22數(shù)值求解舉例數(shù)值求解舉例l 先編寫函數(shù)文件先編寫函數(shù)文件 verderpol.mfunction xprime=verderpol(t,x)global mu;x
16、prime=x(2); mu*(1-x(1)2)*x(2) - x(1);l 再編寫腳本文件再編寫腳本文件 vdpl.m,在命令窗口直接運(yùn)行該文件。,在命令窗口直接運(yùn)行該文件。 clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0); plot(t,x(:,1),r-);23Matlab 求解微分方程小結(jié)求解微分方程小結(jié)q Matlab 函數(shù)函數(shù)u 求解析解(通解或特解),用求解析解(通解或特解),用 dsolveu 求數(shù)值解(特解),用求數(shù)值解(特解),用 ode45、ode23 .q Matlab 編程編程u Euler 折線法折線
17、法u Runga-Kutta 方法方法24*n想:無論是解析解還是數(shù)值解,都不如圖形解直觀明了。在想:無論是解析解還是數(shù)值解,都不如圖形解直觀明了。在MATLAB下用下用 什么命令可以什么命令可以用圖形顯示解析解或數(shù)值解?用圖形顯示解析解或數(shù)值解?n提示:如果微分方程(組)的解析解為提示:如果微分方程(組)的解析解為y=f(x),則可以用,則可以用MATLAB函數(shù)函數(shù)fplot作出函數(shù)作出函數(shù)曲線曲線y=f(x)的圖形,具體用法的圖形,具體用法其中其中fun給出函數(shù)表達(dá)式;給出函數(shù)表達(dá)式;lims=xmin xmax ymin ymax限定坐標(biāo)軸的大小例如限定坐標(biāo)軸的大小例如如果得到微分方程(組)的數(shù)值解(
溫馨提示
- 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請(qǐng)下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請(qǐng)聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會(huì)有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
- 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
- 5. 人人文庫網(wǎng)僅提供信息存儲(chǔ)空間,僅對(duì)用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對(duì)用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對(duì)任何下載內(nèi)容負(fù)責(zé)。
- 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請(qǐng)與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時(shí)也不承擔(dān)用戶因使用這些下載資源對(duì)自己和他人造成任何形式的傷害或損失。
最新文檔
- 江西電力職業(yè)技術(shù)學(xué)院《社會(huì)保障史》2023-2024學(xué)年第一學(xué)期期末試卷
- 吉林體育學(xué)院《解剖和生理(生理)》2023-2024學(xué)年第一學(xué)期期末試卷
- 湖南食品藥品職業(yè)學(xué)院《高等數(shù)理統(tǒng)計(jì)一》2023-2024學(xué)年第一學(xué)期期末試卷
- 湖南大眾傳媒職業(yè)技術(shù)學(xué)院《國際物流與供應(yīng)鏈管理》2023-2024學(xué)年第一學(xué)期期末試卷
- 【物理】《壓強(qiáng)》(教學(xué)設(shè)計(jì))-2024-2025學(xué)年人教版(2024)初中物理八年級(jí)下冊(cè)
- 高考物理總復(fù)習(xí)《力與物體的平衡》專項(xiàng)測(cè)試卷含答案
- 年產(chǎn)5萬臺(tái)模塊化操動(dòng)機(jī)構(gòu)及高壓斷路器生產(chǎn)線項(xiàng)目可行性研究報(bào)告寫作模板-備案審批
- 鄭州軟件職業(yè)技術(shù)學(xué)院《鋼結(jié)構(gòu)設(shè)計(jì)B》2023-2024學(xué)年第一學(xué)期期末試卷
- 浙江電力職業(yè)技術(shù)學(xué)院《光電技術(shù)基礎(chǔ)》2023-2024學(xué)年第一學(xué)期期末試卷
- 長安大學(xué)《個(gè)體防護(hù)》2023-2024學(xué)年第一學(xué)期期末試卷
- 上海車位交易指南(2024版)
- 醫(yī)學(xué)脂質(zhì)的構(gòu)成功能及分析專題課件
- 新疆塔城地區(qū)(2024年-2025年小學(xué)六年級(jí)語文)部編版期末考試(下學(xué)期)試卷及答案
- 2024年9月時(shí)事政治試題帶答案
- 汽車供應(yīng)商審核培訓(xùn)
- 《計(jì)算機(jī)網(wǎng)絡(luò) 》課件第1章
- 1《地球的表面》說課稿-2024-2025學(xué)年科學(xué)五年級(jí)上冊(cè)教科版
- GB/T 44764-2024石油、石化和天然氣工業(yè)腐蝕性石油煉制環(huán)境中抗硫化物應(yīng)力開裂的金屬材料
- 自動(dòng)化招聘筆試試題及答案
- 重慶市主城四區(qū)2025屆高一物理第一學(xué)期期末聯(lián)考試題含解析
- 燃?xì)庑袠I(yè)有限空間作業(yè)安全管理制度
評(píng)論
0/150
提交評(píng)論