實用 matlab教學(xué)資料-27-28常微分方程續(xù)ppt課件

1、適用適用MATLAB剛性問題剛性問題ode solver求 解 器求 解 器Solver算法算法ODE類型類型精度精度適用場合適用場合ode15sGear法法剛性剛性低低中中若因剛性導(dǎo)致若因剛性導(dǎo)致ode45失失效時，可嘗試使用；效時，可嘗試使用；可可求解求解DAEs(微分代數(shù)方微分代數(shù)方程程)ode23s單 步單 步 2 階階Rosenbrock法法剛性剛性低低解決精度要求低、剛性解決精度要求低、剛性問 題 。 計 算 時 間 比問 題 。 計 算 時 間 比ode15s短短ode23t梯形梯形Euler公式公式中等剛性中等剛性低低適用于中等剛性的場合；適用于中等剛性的場合；求解求解DAEs

2、ode23tb隱 式隱 式 R K公式公式剛性剛性低低解決精度要求低、剛性解決精度要求低、剛性問題。精度較低時，計問題。精度較低時，計算時間比算時間比ode15s短短例例4 分別用分別用ode15s和和ode45解如下解如下剛性方程剛性方程1tspan=0,300 2tspan=0,3000 1222121121000(1)(0)2(0)0yyyyyyyy03000 x5、隱式常微分方程、隱式常微分方程( , ,)0(0)(0)f x y yyayb定解條件：函數(shù)和導(dǎo)函數(shù)初值定解條件：函數(shù)和導(dǎo)函數(shù)初值5、隱式常微分方程、隱式常微分方程ode15ix,y=ode15i(odefun,tspan,

3、y0, dydx0,) n定義隱式方程，三個輸入定義隱式方程，三個輸入x,y,y( , , )0f x y yn函數(shù)初值函數(shù)初值n導(dǎo)函數(shù)初導(dǎo)函數(shù)初值值例例5 求隱式常微分方程在求隱式常微分方程在0,12的解的解11221122121212sincos210coscos30(0)0.1(0)0.2(0)0(0)1yyyyyyyyyyyyyy x,y=ode15i(odefun,tspan,y0, dydx0) 6、延遲常微分方程、延遲常微分方程n描畫延遲系統(tǒng)：描畫延遲系統(tǒng)：n 系統(tǒng)隨時間的演化，不僅依賴于系系統(tǒng)隨時間的演化，不僅依賴于系統(tǒng)當(dāng)前的形狀，且依賴于系統(tǒng)過去某一時統(tǒng)當(dāng)前的形狀，且依賴于系

4、統(tǒng)過去某一時辰或假設(shè)干時辰的形狀辰或假設(shè)干時辰的形狀1( ), ( ), (),. ()ny tf t y ty ty t0延遲常數(shù)延遲常數(shù)dde23適用場所：適用場所：解恒定遲滯的延遲微分方程解恒定遲滯的延遲微分方程原理：隱式原理：隱式Runge-Katta算法算法dde23options =ddeset(name1,value1,.)sol=dde23(ddefun,lags,history,tspan,options) lags-延遲常數(shù)，行向量延遲常數(shù)，行向量12l ,.0njags dde23sol=dde23(ddefun,lags,history,tspan,options) f

5、unction dydt = ddefun(x,y,Z) 其中其中z待解函數(shù)的歷史形狀，矩陣待解函數(shù)的歷史形狀，矩陣11211(:,1)();();.()mzy ty tyt12(:, )();();.()jjmjzjy ty tytj = lags(j)j = lags(j)函數(shù)式函數(shù)式m文件定義文件定義dde23sol=dde23(ddefun,lags,history,tspan,options) history-定義歷史形狀，定義歷史形狀，待解函數(shù)待解函數(shù)y(t)在在t t0的值，的值，列向量列向量options =ddeset(name1,value1,.)參數(shù)設(shè)置參數(shù)設(shè)置RelTo

6、l-相對誤差相對誤差A(yù)bsTol-絕對誤差絕對誤差NormControl-范數(shù)誤差控制范數(shù)誤差控制InitialStep-初始步長初始步長MaxStep 最大步長最大步長ddesd適用場所：適用場所：處置遲滯是自變量函數(shù)的延遲微分處置遲滯是自變量函數(shù)的延遲微分方程方程ddesdsol=ddesd(ddefun,lags,history,tspan,options) lags:普通延遲普通延遲,lags=p(t)用匿名函數(shù)或函數(shù)式用匿名函數(shù)或函數(shù)式m文件定義文件定義 例例6 解延遲常微分方程解延遲常微分方程9.642 (3)( )( )1(3)y ty ty ty t0( )0.6ty t，n1

7、 1、邊值問題的描畫、邊值問題的描畫在自變量的兩端給定邊境條件在自變量的兩端給定邊境條件10.2 邊值問題邊值問題n一、邊值問題數(shù)值解法一、邊值問題數(shù)值解法( , ,)( )( )yf x y yaxby ay bnBoundary Value Problemsn有限差分法有限差分法n 將邊值問題包括方程及邊將邊值問題包括方程及邊境條件離散化，求離散點(diǎn)上境條件離散化，求離散點(diǎn)上函數(shù)近似值函數(shù)近似值yi。2、數(shù)值求解方法、數(shù)值求解方法nBoundary Value Problemsbvp4c, bvp5c功能：解常微分方程邊值問題功能：解常微分方程邊值問題算法：有限差分法算法：有限差分法n二、二

8、、MATLAB功能函數(shù)功能函數(shù)nBoundary Value Problemssol=bvp4c(odefun,bcfun,solinit,options)bvp4cnBoundary Value Problemsn定義微分方定義微分方程程n定義邊境條定義邊境條件件n初始猜初始猜測測sol=bvp4c(odefun,bcfun,solinit,options)(1)odefun定義降階后構(gòu)成的常微分方程組定義降階后構(gòu)成的常微分方程組, 同同ode solver必需前往導(dǎo)函數(shù)列向量必需前往導(dǎo)函數(shù)列向量sol=bvp4c(odefun,bcfun,solinit,options)(2) bcfun

9、 定義邊境條件定義邊境條件 res = bcfun(ya,yb) 0101( )( )( )( )a y aa y ab y bb y b0101( )( )0( )( )0a y aa y ab y bb y b res sol=bvp4c(odefun,bcfun,solinit,options)(2) bcfun 定義邊境條件定義邊境條件 res = bcfun(ya,yb)res：邊境條件殘差向量，由：邊境條件殘差向量，由ya,yb 表達(dá)表達(dá)ya, yb 分別為左、右邊境條件分別為左、右邊境條件ya(1), yb(1)表示邊境點(diǎn)函數(shù)值表示邊境點(diǎn)函數(shù)值ya(2), yb(2)表示邊境點(diǎn)一

10、階導(dǎo)函數(shù)值，表示邊境點(diǎn)一階導(dǎo)函數(shù)值， 依次類推。依次類推。(3) solinit：解的初始猜測網(wǎng)格：解的初始猜測網(wǎng)格 經(jīng)過函數(shù)經(jīng)過函數(shù)bvpinit實現(xiàn)實現(xiàn)solinit = bvpinit(x,yinit)其中：其中：x-可由線性間隔向量生成，作為初始網(wǎng)格節(jié)點(diǎn)可由線性間隔向量生成，作為初始網(wǎng)格節(jié)點(diǎn)yinit-初值猜測向量，與微分方程組維數(shù)一致初值猜測向量，與微分方程組維數(shù)一致yinit(1)作為作為y1在一切網(wǎng)格節(jié)點(diǎn)在一切網(wǎng)格節(jié)點(diǎn)x(i)處的初值處的初值yinit(2)作為作為y2在一切網(wǎng)格節(jié)點(diǎn)在一切網(wǎng)格節(jié)點(diǎn)x(i)處的初值處的初值. sol=bvp4c(odefun,bcfun,solin

11、it,options)sol=bvp4c(odefun,bcfun,solinit,options)(4)optionsoptions=bvpset(name1,value1,.)namevaluedescriptionRelToldefault, 1e-3相對誤差相對誤差A(yù)bsTol default, 1e-6絕對誤差絕對誤差誤差控制誤差控制(4)Options數(shù)值解計算信息輸出控制數(shù)值解計算信息輸出控制name value包括信息包括信息Stats on Off(默認(rèn)默認(rèn))網(wǎng)格數(shù)目網(wǎng)格數(shù)目最大殘差最大殘差調(diào)用微分方程的次數(shù)調(diào)用微分方程的次數(shù)調(diào)用邊界條件的次數(shù)調(diào)用邊界條件的次數(shù)options

12、=bvpset(name1,value1,.)nsol.x-實踐計算網(wǎng)格節(jié)點(diǎn)實踐計算網(wǎng)格節(jié)點(diǎn)nsol.y-解向量解向量ynsol.yp-節(jié)點(diǎn)導(dǎo)數(shù)值節(jié)點(diǎn)導(dǎo)數(shù)值ynsol.parameters- 假設(shè)涉及未知參數(shù)求假設(shè)涉及未知參數(shù)求解，那么前往求出的參數(shù)解，那么前往求出的參數(shù)nsol.solver-求解器求解器nsol.stats-計算統(tǒng)計數(shù)據(jù)計算統(tǒng)計數(shù)據(jù)(5)sol: 解的構(gòu)造數(shù)組解的構(gòu)造數(shù)組sol=bvp4c(odefun,bcfun,solinit,options)其中其中x-輸出網(wǎng)格節(jié)點(diǎn)輸出網(wǎng)格節(jié)點(diǎn)idx-正整數(shù)正整數(shù),表示要輸出數(shù)值解的函數(shù)序號表示要輸出數(shù)值解的函數(shù)序號,例例如如 2表示

13、輸出表示輸出y2的數(shù)值解。假設(shè)不寫的數(shù)值解。假設(shè)不寫,輸出一輸出一切函數(shù)數(shù)值解切函數(shù)數(shù)值解sx-數(shù)值解數(shù)值解 sx(1,:)是是y1的數(shù)值解的數(shù)值解 sx(:,1)是一切函數(shù)在第一個節(jié)點(diǎn)的數(shù)值解是一切函數(shù)在第一個節(jié)點(diǎn)的數(shù)值解(6)deval:獲取數(shù)值解獲取數(shù)值解sx = deval(sol, x, idx)sol=bvp4c(odefun,bcfun,solinit,options)options =bvpset(name1,value1,.)%參數(shù)設(shè)置參數(shù)設(shè)置solinit = bvpinit(x,yinit)%初始網(wǎng)格和初值的設(shè)置初始網(wǎng)格和初值的設(shè)置sol= bvp4c(odefun,bc

14、fun,solinit,options)sx = deval(sol, x, idx) %獲取數(shù)值解獲取數(shù)值解bvp4c的普通格式的普通格式nBoundary Value Problems節(jié)節(jié) 點(diǎn)點(diǎn)含含 義義設(shè)設(shè) 置置workspace中的中的變量名變量名初始網(wǎng)格初始網(wǎng)格節(jié)點(diǎn)節(jié)點(diǎn)用于計算網(wǎng)格用于計算網(wǎng)格初值初值solinit = bvpinit(x,yinit)結(jié) 構(gòu) 數(shù) 組結(jié) 構(gòu) 數(shù) 組solint中的中的solint.x實際計算實際計算網(wǎng)格節(jié)點(diǎn)網(wǎng)格節(jié)點(diǎn)MATLAB計算計算時采用的時采用的節(jié)點(diǎn)節(jié)點(diǎn)軟件根據(jù)精度要求自行軟件根據(jù)精度要求自行設(shè)置，用戶無法修改。設(shè)置，用戶無法修改。結(jié)構(gòu)數(shù)組結(jié)構(gòu)數(shù)組Sol中的中的Sol.x輸出網(wǎng)格輸出網(wǎng)格節(jié)點(diǎn)節(jié)點(diǎn)根據(jù)用戶需要，根據(jù)用戶需要，對輸出數(shù)據(jù)的對輸出數(shù)據(jù)的疏密進(jìn)行設(shè)置疏密進(jìn)行設(shè)置sx = deval(sol,x,idx