重慶工商大學(xué)數(shù)學(xué)模型與數(shù)學(xué)實(shí)驗(yàn)課件第08講 常微分方程數(shù)值解

1、大學(xué)數(shù)學(xué)實(shí)驗(yàn)Experiments in Mathematics實(shí)驗(yàn)4 常微分方程數(shù)值解數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)- 李煥榮為什么要學(xué)習(xí)微分方程數(shù)值解 微分方程是研究函數(shù)變化規(guī)律的重要工具，有著廣泛的應(yīng)用。如: 物體的運(yùn)動(dòng), 電路的電壓, 人口增長的預(yù)測(cè) 許多微分方程沒有解析解，數(shù)值解法是求解的重要手段，如實(shí)驗(yàn)4的基本內(nèi)容3. 實(shí)際問題用微分方程建模，并求解2. 龍格-庫塔方法的MATLAB實(shí)現(xiàn)*4. 數(shù)值算法的收斂性、穩(wěn)定性與剛性方程兩個(gè)最常用的數(shù)值算法: 歐拉（Euler）方法 龍格-庫塔（Runge-Kutta）方法實(shí)例1 海上緝私 海防某部緝私艇上的雷達(dá)發(fā)現(xiàn)正東方向c海里處有一艘走私船正以速度

2、a向正北方向行駛，緝私艇立即以最大速度b(a)前往攔截。如果用雷達(dá)進(jìn)行跟蹤時(shí)，可保持緝私艇的速度方向始終指向走私船。 建立任意時(shí)刻緝私艇位置及 航線的數(shù)學(xué)模型,并求解; 求出緝私艇追上走私船的時(shí)間。 a北bc艇船實(shí)例1 海上緝私 建立坐標(biāo)系如圖: t=0 艇在(0, 0), 船在(c, 0); 船速a, 艇速b 時(shí)刻 t 艇位于P(x, y), 船到達(dá) Q(c, at)模型:0yxcR(c,y )Q(c,at)P(x,y)b由方程無法得到x(t), y(t)的解析解需要用數(shù)值解法求解 “常微分方程初值問題數(shù)值解”的提法不求解析解 y = y(x) (無解析解或求解困難)而在一系列離散點(diǎn)通常取等

3、步長hy1y2ynyP0 x0 x1x2x3 xy=y(x)y0P1P2P3歐 拉 方 法基本思路x取不同點(diǎn)各種歐拉公式向前歐拉公式顯式公式P3歐 拉 方 法向后歐拉公式 隱式公式y(tǒng)P0 x0 x1x2x3 xy0y=y(x)P1P2向前歐拉公式向后歐拉公式二者平均得到梯形公式仍為隱式公式，需迭代求解將梯形公式的迭代過程簡化為兩步預(yù)測(cè)校正改進(jìn)歐拉公式龍格-庫塔方法 向前，向后歐拉公式： 用xn, xn+1內(nèi)個(gè)點(diǎn)的導(dǎo) 數(shù)代替 f(x, y(x) 梯形公式，改進(jìn)歐拉公式：用xn, xn+1內(nèi)個(gè)點(diǎn)導(dǎo)數(shù)的平均值代替 f(x, y(x)龍格-庫塔方法的基本思想在xn, xn+1內(nèi)多取幾個(gè)點(diǎn)，將它們的導(dǎo)數(shù)

4、加權(quán)平均代替 f(x, y(x)，設(shè)法構(gòu)造出精度更高的計(jì)算公式。 常用的(經(jīng)典)龍格庫塔公式不足:收斂速度較慢L級(jí)?階龍格-庫塔方法的一般形式使精度盡量高4級(jí)4階微分方程組和高階方程初值問題的數(shù)值解 歐拉方法和龍格-庫塔方法可直接推廣到微分方程組向前歐拉公式 高階方程需要先降階為一階微分方程組 龍格庫塔方法的 MATLAB 實(shí)現(xiàn) t,x=ode23(f, ts, x0) 3級(jí)2階龍格-庫塔公式 t,x=ode45(f, ts, x0) 5級(jí)4階龍格-庫塔公式 f是待解方程寫成的函數(shù)m文件： function dx=f(t,x)dx=f1; f2; fn;ts = t0,t1, ,tf輸出指定時(shí)

5、刻 t0,t1, ,tf 的函數(shù)值ts = t0:k:tf輸出 t0,tf 內(nèi)等分點(diǎn)處的函數(shù)值x0為函數(shù)初值(n維)輸出t=ts, x為相應(yīng)函數(shù)值(n維) 缺省精度(相對(duì)誤差10-3，絕對(duì)誤差10-6), 計(jì)算步長按精度要求自動(dòng)調(diào)整.實(shí)例1海上緝私(續(xù)) 模型的數(shù)值解0yxc(x(t),y(t)ab設(shè): 船速a=20 (海里/小時(shí)) 艇速b=40 (海里/小時(shí)) 距離c=15 (海里) 求: 緝私艇的位置x(t),y(t) 緝私艇的航線y(x) 實(shí)例1海上緝私(續(xù)) 模型的數(shù)值解 a=20, b=40, c=15走私船的位置x1(t)= c=15y1(t)=at=20t t=0.5時(shí)緝私艇追上

6、走私船 緝私艇的航線y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.94450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.05.06.07.08.09.010.0 剛性現(xiàn)象與剛性方程振動(dòng)系統(tǒng)或電路系統(tǒng)的數(shù)學(xué)模型 現(xiàn)象 k=2000.5, r=1000, a=1, b=-1999.5, f(t)=1 瞬態(tài)解與穩(wěn)態(tài)解 e-2000t 快瞬態(tài)解 e-t/2 慢瞬態(tài)解 計(jì)算到t=0.005時(shí)已衰減到4.510-5 計(jì)算到t=20時(shí)才衰減到4.510-5 精度達(dá)到10-4需算到t=20(由慢瞬態(tài)解=1/2決定)選取步長h由快瞬態(tài)解=2000決定h2.785/2000=0.0014 龍格-庫塔公式t=20需14286步 快、慢瞬態(tài)解的特征根相差懸殊 剛性現(xiàn)象(Stiff) 求穩(wěn)態(tài)解 剛性現(xiàn)象與剛性方程剛性方程的MATLAB求解 ode23, ode45 解剛性方程的困難步長自動(dòng)變小 計(jì)算時(shí)間很長 求解剛性方程的命令: ode23s, ode15s 等 (用法相同)例特征根1