第七章微積分的數(shù)值計(jì)算

1、7.1 數(shù)值微分7.2 數(shù)值積分7.3 常微分方程的數(shù)值解法 7.1 數(shù)值微分?jǐn)?shù)值微分 實(shí)際問題常需計(jì)算函數(shù)的導(dǎo)數(shù)或積分值。但很多情況下，函數(shù)關(guān)系難以準(zhǔn)確表示；即使能使用解析式準(zhǔn)確表示，表示式卻很復(fù)雜，不能用于實(shí)際計(jì)算。本章介紹數(shù)值計(jì)算導(dǎo)數(shù)或積分的實(shí)用方法。7.1.1 差分和差商差分和差商 根據(jù)導(dǎo)數(shù)的定義其中，x和y分別稱為自變量x和因變量y的增量，也稱之為差分?？梢杂貌罘值纳?作微商（導(dǎo)數(shù)）的近似。數(shù)值微分就是用函數(shù)值的線性組合近似函數(shù)在某點(diǎn)的導(dǎo)數(shù)值。自變量x的步長一般取定值。 首先在xi處對(duì)函數(shù)進(jìn)行泰勒展開，1001limlimiiixxiiyyyyxxx yx234( 4 )12 !3

2、!4 !iiiiiixxxyyxyyyy 根據(jù)不同的組合方式可以得到精度不同的差分公式。以函數(shù)的一階導(dǎo)數(shù)為例 ：微分公式微分公式表達(dá)式表達(dá)式截?cái)嗾`差截?cái)嗾`差兩點(diǎn)前向o(x)兩點(diǎn)后向o(x)三點(diǎn)中心o(x2)三點(diǎn)前向o(x2)三點(diǎn)后向o(x2)五點(diǎn)中心o(x4)1iiiyyyx 1iiiyyyx 112iiiyyyx 12432iiiiyyyyx12432iiiiyyyyx 11228812iiiiiyyyyyx精度為o(x2)的高階中心差分算法11221123(4)211242222464iiiiiiiiiiiiiiiyyyyxyyyyyxyyyyyyx精度為o(x4)的高階中心差分算法211

3、2211223211233(4)3211234881216301612813138812395639126iiiiiiiiiiiiiiiiiiiiiiiiiiyyyyyxyyyyyyxyyyyyyyxyyyyyyyyx7.1.2 數(shù)值微分的實(shí)現(xiàn)數(shù)值微分的實(shí)現(xiàn) 在matlab中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)diff和梯度函數(shù)gradient。diff調(diào)用格式為：調(diào)用格式為： dy=diff(y)：計(jì)算向量y的向前差分，并把結(jié)果賦值給向量dy dy(i)=y(i+1)-y(i)，i=1,2,n-1。注意向量dy元素個(gè)數(shù)比y少dy=diff(y,n)：計(jì)算向量y的n階向前差分

4、。例如， diff(y,2)=diff(diff(y)=dx(i+1)-dx(i)= y(i+2)-2y(i+1)+y(i) ， i=1,2 n-2。dx=diff(a,n,dim)：計(jì)算矩陣a的n階差分，dim=1時(shí)(缺省狀態(tài))，按列計(jì)算差分；dim=2，按行計(jì)算差分。 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 b=diff(a)b = 0 1 2 3 0 1 3 6 0 1 4 10 c=diff(a,2)c = 0 0 1 3 0 0 1 4 d=diff(a,1,1)d = 0 1 2 3 0 1 3 6 0 1 4 10 e=

5、diff(a,1,2)e = 0 0 0 1 1 1 2 3 4 3 6 107.1.3 近似梯度函數(shù)近似梯度函數(shù)gradient的調(diào)用格式為的調(diào)用格式為 fx,fy=gradient(f,hx,hy) 求矩陣f， 求其x(行)方向的數(shù)值梯度fx和y(列)方向的數(shù)值梯度fy，x方向步長全為hx，y方向步長全為hy。 fx相當(dāng)于偏導(dǎo)數(shù)f/x ,fy相當(dāng)于偏導(dǎo)數(shù)f/y。 fx,fy=gradient(f,h) 求矩陣f， 求其x(行)方向的數(shù)值梯度fx和y(列)方向的數(shù)值梯度fy，各個(gè)方向步長全為h。 fx=gradient(f) 如果f是向量，直接求其數(shù)值梯度；如果f是矩陣，求其x(行)方向的數(shù)

6、值梯度，步長為1 。 fx,fy=gradient(f) 求矩陣f， 求其x(行)方向的數(shù)值梯度fx和y(列)方向的數(shù)值梯度fy，步長為1ffffijkxyz x=1 3 5 2 4;y= gradient(x)y = 2.0000 2.0000 -0.5000 -0.5000 2.0000 y=gradient(x,2)y= 1.0000 1.0000 -0.2500 -0.2500 1.0000即兩兩邊用前向和后向差分，中間用中邊用前向和后向差分，中間用中心差分心差分 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 fx,fy=gradi

7、ent(a)fx = 0 0 0 0 1.0000 1.0000 1.0000 1.0000 2.0000 2.5000 3.5000 4.0000 3.0000 4.5000 8.0000 10.0000fy = 0 1.0000 2.0000 3.0000 0 1.0000 2.5000 4.5000 0 1.0000 3.5000 8.0000 0 1.0000 4.0000 10.0000211315342234545,222xxyhxxxxxxxyyyhxhxhxxxyhx，7.1.4 拉普拉斯算子拉普拉斯算子4*del2 由于內(nèi)部算法的原因： u=4*del2(v,h)，對(duì)1維向量

8、v以步長h求拉普拉斯算子時(shí)，返回一相同維數(shù)的向量u，且默認(rèn)的步長為1。 u=4*del2(v,h1,h2)，對(duì)矩陣v，橫向(x方向)以步長h1，縱向(y方向)以步長 h2計(jì)算拉普拉斯算子。2222211( ,)( ,)( ,)44vx yvx yuvx yxy14321211212321(452)1(2)1( 254)iiiinnnnnuvvvvhuvvvhuvvvvh 4*del2(u)ans = 4 4 4 4 4 4 4 4 4 4 4 4x,y= meshgrid(1:4,1:3); u = x.*x+y.*yu = 2 5 10 17 5 8 13 20 10 13 18 257.2

9、 數(shù)值積分?jǐn)?shù)值積分7.2.1 數(shù)值積分基本原理數(shù)值積分基本原理 我 們 知 道 ， 定 積 分 是 求 和 式 的 極 限 ，即 。它的幾何意義是曲邊梯形的面積。從定義可知，定積分的基本分析方法是四步，即分割、近似、求和、取極限。分割就是把總量（整塊曲邊梯形面積）分成若干分量（小曲邊梯形面積）；近似就是在每個(gè)分量中用容易計(jì)算的量去代表；求和就是把分量加起來得到總近似值；最后取極限就得到積分精確值。 1( )dlim()nbkkankfxxfx把區(qū)間a，b分割成n等分，像這樣取定步長算積分的方法，稱為定步長積分法法。常見的低階求積分公式常見的低階求積分公式 復(fù)化矩形公式 復(fù)化的梯形公式 復(fù)化的辛

10、普森（simpson）公式 10( )bnkkaf x dxs1/21()4 ()()6kkkkhsf xf xf x()kksfxh1()()2kkkhsfxfx, (0,1,1)kbaxakh hknn辛普森公式的幾何意義辛普森公式的幾何意義 7.2.2 變步長積分法變步長積分法 計(jì)算積分，可以采取逐步縮小步長h的辦法。即先任取步長h進(jìn)行計(jì)算，然后取較小步長 h 進(jìn)行計(jì)算，如果兩次計(jì)算結(jié)果相差較大，則取更小步長進(jìn)行計(jì)算，如此下去，直到相鄰兩次計(jì)算結(jié)果相差不大為止，取最小步長算出的結(jié)果作為積分值。這種方法稱為變步長積分法。 利用兩種步長計(jì)算積分時(shí)，通常取h=h/2 。而每次改變步長后，只需計(jì)

11、算新增節(jié)點(diǎn)處的函數(shù)值，將它們的和乘新步長 。 matlab常用的數(shù)值積分函數(shù)參數(shù)名參數(shù)名功能描述功能描述quad一元函數(shù)的數(shù)值積分，采用變步長辛普森simpson法(低階)quadl一元函數(shù)的數(shù)值積分，采用變步長洛巴托l(wèi)obatto法(高階)qualv一元函數(shù)的矢量數(shù)值積分dbquad二重積分（默認(rèn)值為simpson法）triplequad三重積分（默認(rèn)值為simpson法） 7.2.3 一元函數(shù)的數(shù)值積分舉例一元函數(shù)的數(shù)值積分舉例求定積分1. 建立匿名函數(shù)建立匿名函數(shù) f=(x)x./(x.2+1);2. 用用sum函數(shù)實(shí)現(xiàn)復(fù)化矩形法求積函數(shù)實(shí)現(xiàn)復(fù)化矩形法求積1122001lo g (1)0

12、 .3 4 6 612exd xxxx=0:0.001:1; %步長提高到0.001 y=f(x); sum(y)*0.001ans = 0.3468 x=0:0.01:1; %步長0.01 y=f(x); sum(y)*0.01ans = 0.3491 3. 用用trapz函數(shù)實(shí)現(xiàn)復(fù)化梯形法求積函數(shù)實(shí)現(xiàn)復(fù)化梯形法求積x=0:0.001:1;y=f(x);y=trapz(y)*0.001y = 0.3466x=0:0.01:1;y=f(x);y=trapz(y)*0.01y = 0.34664. 用用matlab函數(shù)求定積分函數(shù)求定積分用quad實(shí)現(xiàn)變步長辛普森法求定積分。調(diào)用格式為 quad

13、(fun,a,b,tol)其中fun為積分函數(shù)，a,b為積分區(qū)間，tol為積分的誤差閾值，默認(rèn)值為1e-6quad(f,0,1)ans = 0.3466用quadl實(shí)現(xiàn)變步長洛巴托求定積分。調(diào)用格式為 quadl(fun,a,b,tol)其中fun為積分函數(shù)，a,b為積分區(qū)間，tol為積分的誤差閾值，默認(rèn)值為1e-6 quadl(f,0,1)ans = 0.34667.2.4 矢量積分矢量積分相當(dāng)于多個(gè)一元定積分。例如求y=(x,n)1./(sqrt(2*pi).*(1:n).*exp(-x.2./(2*(1:n).2); %歸一化高斯函數(shù)quadv(x)y(x,5),-1,1)ans = 0

14、.6827 0.3829 0.2611 0.1974 0.1585矢量積分的結(jié)果是一個(gè)向量，每一個(gè)元素為一個(gè)一元函數(shù)定積分的值。12211exp(),1,2,3,4,522xdx nn7.2.5 二重定積分的數(shù)值求解二重定積分的數(shù)值求解 使用matlab提供的dblquad函數(shù)就可以直接求出上述二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為： dblquad(f,a,b,c,d,tol,method)該函數(shù)求f(x,y)在a,bc,d區(qū)域上的二重定積分。參數(shù)tol，可以采用method=quadl的方法使函數(shù)用高階的洛巴托法計(jì)算定積分。例如計(jì)算 f=(x,y)exp(-x.2-y.2)/pi; %歸一

15、化高斯函數(shù)dblquad(f,-1,1,-1,1,1e-6,quadl)ans = 0.7101三重積分函數(shù)triplequad用法與二重積分類似1122111exp()xyd xdy 7.3 常微分方程的數(shù)值解法常微分方程的數(shù)值解法7.3.1 常微分方程初值問題的數(shù)值解法常微分方程初值問題的數(shù)值解法一般形式為所謂的數(shù)值解法就是求解y(x)在 區(qū)間的近似值yn的方法，yn(n=1,2,n)稱為常微分方程的數(shù)值解。自變量x的步長一般為定值h。0( ,)( )dyfx yaxbdxy ay01na xxxb7.3.2 常見的數(shù)值方法常見的數(shù)值方法向前歐拉公式向后歐拉公式梯形公式改進(jìn)的歐拉公式10(

16、,)0,1( )nnnnyyhfxynnyy a1110(,)0,1( )nnnnyyhfxynnyy a1110(,)(,)0,12( )nnnnnnhyyf xyf xynnyy a1(,)(,)1()2pnnnqnnpnpqyyhfxyyyhfxh yyyy 7.3.3 龍格庫塔法簡(jiǎn)介龍格庫塔法簡(jiǎn)介 基本思想就是利用在某些點(diǎn)處的值的線性組合構(gòu)造公式，使其按泰勒展開后與初值問題的解的泰勒展開式比較，有盡可能多的相同項(xiàng)，從而保證算式有較高的精度。常用的四階經(jīng)典的龍格-庫塔公式112341213243(22)6(,)(,)22(,)22(,)22nnnnnnnnnnhyykkkkkf x yh

17、hkf xykhhkf xykhhkf xyk7.3.4 龍格庫塔法的實(shí)現(xiàn)龍格庫塔法的實(shí)現(xiàn)求解器求解器ode類型類型特點(diǎn)特點(diǎn)精度精度說明說明ode45非剛性一步算法（只需前一步的結(jié)果），4,5階 runge-kutta方法。中大部分場(chǎng)合的首選算法ode23非剛性一步算法，2,3階runge-kutta方法。低使用于精度較低的情形ode113非剛性多步法（需要前幾布的結(jié)果），adams算法。低高計(jì)算時(shí)間比ode45短ode23t適度剛性 采用梯形算法適度剛性情形ode15s剛性多步法，gears反向數(shù)值積分。低中若ode45失效時(shí)，可嘗試使用ode23s剛性一步法，2階rosebrock算法精度

18、。低當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短 對(duì)于一個(gè)常微分方程組，如果其解相差十分懸殊，就稱之為剛性方程組剛性方程組。對(duì)于剛性方程組，為了保持解法的穩(wěn)定，步長選取十分困難，有些解法不再適用。ode函數(shù)調(diào)用格式：t,y = ode ij (odefun,tspan,y0)t,y = ode ij (odefun,tspan,y0,options)t,y = ode ij (odefun,tspan,y0,options,p1,p2,.)t,y,te,ye,ie=odeij(odefun,tspan,y0,options,p1,p2,.)odefun為顯式常微分方程中的 f(xn,yn)，tspan為求解區(qū)間，y0為初始條件。7.3.5 求解多變量一階常微分方程組。求解多變量一階常微分方程組。例：求解下面的混沌理論的洛侖茲方程。dx8= -y + yzdt3dy= 10y + 10zdtdz= - xy + 35y- zdt編寫lorfun.m如下：function ydot=lorfun(t,y)ydot=-8/3*y(1)+y(2)*y