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

第七章微積分的數(shù)值計(jì)算第一頁(yè)，共三十三頁(yè)，2022年，8月28日

7.1數(shù)值微分實(shí)際問(wèn)題常需計(jì)算函數(shù)的導(dǎo)數(shù)或積分值。但很多情況下，函數(shù)關(guān)系難以準(zhǔn)確表示；即使能使用解析式準(zhǔn)確表示，表示式卻很復(fù)雜，不能用于實(shí)際計(jì)算。本章介紹數(shù)值計(jì)算導(dǎo)數(shù)或積分的實(shí)用方法。第二頁(yè)，共三十三頁(yè)，2022年，8月28日7.1.1差分和差商

根據(jù)導(dǎo)數(shù)的定義其中，?x和?y分別稱為自變量x和因變量y的增量，也稱之為差分。可以用差分的商作微商（導(dǎo)數(shù)）的近似。數(shù)值微分就是用函數(shù)值的線性組合近似函數(shù)在某點(diǎn)的導(dǎo)數(shù)值。自變量x的步長(zhǎng)一般取定值。首先在xi處對(duì)函數(shù)進(jìn)行泰勒展開(kāi)，第三頁(yè)，共三十三頁(yè)，2022年，8月28日根據(jù)不同的組合方式可以得到精度不同的差分公式。以函數(shù)的一階導(dǎo)數(shù)為例：微分公式表達(dá)式截?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)………………第四頁(yè)，共三十三頁(yè)，2022年，8月28日精度為O(?X2)的高階中心差分算法精度為O(?X4)的高階中心差分算法第五頁(yè)，共三十三頁(yè)，2022年，8月28日7.1.2數(shù)值微分的實(shí)現(xiàn)在MATLAB中，沒(méi)有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)diff和梯度函數(shù)gradient。diff調(diào)用格式為：

Dy=diff(Y)：計(jì)算向量Y的向前差分，并把結(jié)果賦值給向量DyDy(i)=Y(i+1)-Y(i)，i=1,2,…,n-1。注意向量Dy元素個(gè)數(shù)比Y少Dy=diff(Y,n)：計(jì)算向量Y的n階向前差分。例如，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ì)算差分。第六頁(yè)，共三十三頁(yè)，2022年，8月28日>>A=pascal(4)A=1111123413610141020>>B=diff(A)B=0123013601410>>C=diff(A,2)C=00130014>>D=diff(A,1,1)D=0123013601410>>E=diff(A,1,2)E=0001112343610第七頁(yè)，共三十三頁(yè)，2022年，8月28日7.1.3近似梯度函數(shù)gradient的調(diào)用格式為[Fx,Fy]=gradient(F,hx,hy)求矩陣F，求其x(行)方向的數(shù)值梯度Fx和y(列)方向的數(shù)值梯度Fy，x方向步長(zhǎng)全為hx，y方向步長(zhǎng)全為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è)方向步長(zhǎng)全為h。Fx=gradient(F)如果F是向量，直接求其數(shù)值梯度；如果F是矩陣，求其x(行)方向的數(shù)值梯度，步長(zhǎng)為1。[Fx,Fy]=gradient(F)求矩陣F，求其x(行)方向的數(shù)值梯度Fx和y(列)方向的數(shù)值梯度Fy，步長(zhǎng)為1第八頁(yè)，共三十三頁(yè)，2022年，8月28日>>X=[13524];>>Y=gradient(X)Y=2.00002.0000-0.5000-0.50002.0000>>Y=gradient(X,2)Y=1.00001.0000-0.2500-0.25001.0000即兩邊用前向和后向差分，中間用中心差分>>A=pascal(4)A=1111123413610141020>>[Fx,Fy]=gradient(A)Fx=00001.00001.00001.00001.00002.00002.50003.50004.00003.00004.50008.000010.0000Fy=01.00002.00003.000001.00002.50004.500001.00003.50008.000001.00004.000010.0000第九頁(yè)，共三十三頁(yè)，2022年，8月28日7.1.4拉普拉斯算子4*del2由于內(nèi)部算法的原因：U=4*del2(v,h)，對(duì)1維向量v以步長(zhǎng)h求拉普拉斯算子時(shí)，返回一相同維數(shù)的向量U，且默認(rèn)的步長(zhǎng)為1。

U=4*del2(v,h1,h2)，對(duì)矩陣v，橫向(x方向)以步長(zhǎng)h1，縱向(y方向)以步長(zhǎng)h2計(jì)算拉普拉斯算子。第十頁(yè)，共三十三頁(yè)，2022年，8月28日>>4*del2(U)ans=444444444444>>[x,y]=meshgrid(1:4,1:3);>>U=x.*x+y.*yU=25101758132010131825第十一頁(yè)，共三十三頁(yè)，2022年，8月28日7.2數(shù)值積分7.2.1數(shù)值積分基本原理

我們知道，定積分是求和式的極限，即。它的幾何意義是曲邊梯形的面積。從定義可知，定積分的基本分析方法是四步，即分割、近似、求和、取極限。分割就是把總量（整塊曲邊梯形面積）分成若干分量（小曲邊梯形面積）；近似就是在每個(gè)分量中用容易計(jì)算的量去代表；求和就是把分量加起來(lái)得到總近似值；最后取極限就得到積分精確值。

第十二頁(yè)，共三十三頁(yè)，2022年，8月28日把區(qū)間［a，b］分割成n等分，像這樣取定步長(zhǎng)算積分的方法，稱為定步長(zhǎng)積分法法。常見(jiàn)的低階求積分公式復(fù)化矩形公式復(fù)化的梯形公式復(fù)化的辛普森（Simpson）公式

第十三頁(yè)，共三十三頁(yè)，2022年，8月28日辛普森公式的幾何意義

第十四頁(yè)，共三十三頁(yè)，2022年，8月28日7.2.2變步長(zhǎng)積分法計(jì)算積分，可以采取逐步縮小步長(zhǎng)h的辦法。即先任取步長(zhǎng)h進(jìn)行計(jì)算，然后取較小步長(zhǎng)h’

進(jìn)行計(jì)算，如果兩次計(jì)算結(jié)果相差較大，則取更小步長(zhǎng)進(jìn)行計(jì)算，如此下去，直到相鄰兩次計(jì)算結(jié)果相差不大為止，取最小步長(zhǎng)算出的結(jié)果作為積分值。這種方法稱為變步長(zhǎng)積分法。

利用兩種步長(zhǎng)計(jì)算積分時(shí)，通常取h’=h/2。而每次改變步長(zhǎng)后，只需計(jì)算新增節(jié)點(diǎn)處的函數(shù)值，將它們的和乘新步長(zhǎng)。

第十五頁(yè)，共三十三頁(yè)，2022年，8月28日MATLAB常用的數(shù)值積分函數(shù)參數(shù)名功能描述quad一元函數(shù)的數(shù)值積分，采用變步長(zhǎng)辛普森Simpson法(低階)quadl一元函數(shù)的數(shù)值積分，采用變步長(zhǎng)洛巴托Lobatto法(高階)qualv一元函數(shù)的矢量數(shù)值積分dbquad二重積分（默認(rèn)值為Simpson法）triplequad三重積分（默認(rèn)值為Simpson法）第十六頁(yè)，共三十三頁(yè)，2022年，8月28日

7.2.3一元函數(shù)的數(shù)值積分舉例

求定積分

1.建立匿名函數(shù)

f=@(x)x./(x.^2+1);2.用sum函數(shù)實(shí)現(xiàn)復(fù)化矩形法求積

x=[0:0.001:1];%步長(zhǎng)提高到0.001>>y=f(x);>>sum(y)*0.001ans=0.3468>>x=[0:0.01:1];%步長(zhǎng)0.01

>>y=f(x);

>>sum(y)*0.01

ans=0.3491

第十七頁(yè)，共三十三頁(yè)，2022年，8月28日3.用trapz函數(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.3466第十八頁(yè)，共三十三頁(yè)，2022年，8月28日4.用MATLAB函數(shù)求定積分用quad實(shí)現(xiàn)變步長(zhǎng)辛普森法求定積分。調(diào)用格式為quad(fun,a,b,tol)其中fun為積分函數(shù)，[a,b]為積分區(qū)間，tol為積分的誤差閾值，默認(rèn)值為1e-6quad(f,0,1)ans=0.3466用quadl實(shí)現(xiàn)變步長(zhǎng)洛巴托求定積分。調(diào)用格式為quadl(fun,a,b,tol)其中fun為積分函數(shù)，[a,b]為積分區(qū)間，tol為積分的誤差閾值，默認(rèn)值為1e-6

quadl(f,0,1)ans=0.3466第十九頁(yè)，共三十三頁(yè)，2022年，8月28日7.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.68270.38290.26110.19740.1585矢量積分的結(jié)果是一個(gè)向量，每一個(gè)元素為一個(gè)一元函數(shù)定積分的值。第二十頁(yè)，共三十三頁(yè)，2022年，8月28日7.2.5二重定積分的數(shù)值求解使用MATLAB提供的dblquad函數(shù)就可以直接求出上述二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為：dblquad(f,a,b,c,d,tol,method)該函數(shù)求f(x,y)在[a,b]×[c,d]區(qū)域上的二重定積分。參數(shù)tol，可以采用method=@quadl的方法使函數(shù)用高階的洛巴托法計(jì)算定積分。第二十一頁(yè)，共三十三頁(yè)，2022年，8月28日例如計(jì)算f=@(x,y)exp(-x.^2-y.^2)/pi;%歸一化高斯函數(shù)dblquad(f,-1,1,-1,1,1e-6,@quadl)ans=0.7101三重積分函數(shù)triplequad用法與二重積分類似第二十二頁(yè)，共三十三頁(yè)，2022年，8月28日7.3常微分方程的數(shù)值解法7.3.1常微分方程初值問(wèn)題的數(shù)值解法一般形式為所謂的數(shù)值解法就是求解y(x)在區(qū)間的近似值yn的方法，yn(n=1,2,……N)稱為常微分方程的數(shù)值解。自變量x的步長(zhǎng)一般為定值h。第二十三頁(yè)，共三十三頁(yè)，2022年，8月28日7.3.2常見(jiàn)的數(shù)值方法向前歐拉公式向后歐拉公式梯形公式改進(jìn)的歐拉公式第二十四頁(yè)，共三十三頁(yè)，2022年，8月28日7.3.3龍格－庫(kù)塔法簡(jiǎn)介基本思想就是利用在某些點(diǎn)處的值的線性組合構(gòu)造公式，使其按泰勒展開(kāi)后與初值問(wèn)題的解的泰勒展開(kāi)式比較，有盡可能多的相同項(xiàng)，從而保證算式有較高的精度。常用的四階經(jīng)典的龍格-庫(kù)塔公式第二十五頁(yè)，共三十三頁(yè)，2022年，8月28日7.3.4龍格－庫(kù)塔法的實(shí)現(xiàn)求解器ODE類型特點(diǎn)精度說(shuō)明ode45非剛性一步算法（只需前一步的結(jié)果），4,5階Runge-Kutta方法。中大部分場(chǎng)合的首選算法ode23非剛性一步算法，2,3階Runge-Kutta方法。低使用于精度較低的情形ode113非剛性多步法（需要前幾布的結(jié)果），Adams算法。低~高計(jì)算時(shí)間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法，Gear’s反向數(shù)值積分。低~中若ode45失效時(shí)，可嘗試使用ode23s剛性一步法，2階Rosebrock算法精度。低當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短第二十六頁(yè)，共三十三頁(yè)，2022年，8月28日對(duì)于一個(gè)常微分方程組，如果其解相差十分懸殊，就稱之為剛性方程組。對(duì)于剛性方程組，為了保持解法的穩(wěn)定，步長(zhǎng)選取十分困難，有些解法不再適用。ode函數(shù)調(diào)用格式：[t,y]=odeij(@odefun,tspan,y0)[t,y]=odeij(@odefun,tspan,y0,options)[t,y]=odeij(@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為初始條件。