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

1、7.1 數(shù)值微分7.2 數(shù)值積分7.3 常微分方程的數(shù)值解法 7.1 數(shù)值微分?jǐn)?shù)值微分 實(shí)際問題常需計(jì)算函數(shù)的導(dǎo)數(shù)或積分值。但很實(shí)際問題常需計(jì)算函數(shù)的導(dǎo)數(shù)或積分值。但很多情況下，函數(shù)關(guān)系難以準(zhǔn)確表示；即使能使用解多情況下，函數(shù)關(guān)系難以準(zhǔn)確表示；即使能使用解析式準(zhǔn)確表示，表示式卻很復(fù)雜，不能用于實(shí)際計(jì)析式準(zhǔn)確表示，表示式卻很復(fù)雜，不能用于實(shí)際計(jì)算。本章介紹數(shù)值計(jì)算導(dǎo)數(shù)或積分的實(shí)用方法。算。本章介紹數(shù)值計(jì)算導(dǎo)數(shù)或積分的實(shí)用方法。7.1.1 差分和差商差分和差商 根據(jù)導(dǎo)數(shù)的定義根據(jù)導(dǎo)數(shù)的定義其中，其中，x和和y分別稱為自變量分別稱為自變量x和因變量和因變量y的的增量，也稱之為差分?？梢杂貌罘值纳淘隽?/p>

2、，也稱之為差分。可以用差分的商 作微商導(dǎo)數(shù)的近似。數(shù)值微分就是用作微商導(dǎo)數(shù)的近似。數(shù)值微分就是用函數(shù)值的線性組合近似函數(shù)在某點(diǎn)的導(dǎo)數(shù)函數(shù)值的線性組合近似函數(shù)在某點(diǎn)的導(dǎo)數(shù)值。自變量值。自變量x的步長一般取定值。的步長一般取定值。 首先在首先在xi處對函數(shù)進(jìn)行泰勒展開，處對函數(shù)進(jìn)行泰勒展開，1001limlimiiixxiiyyyyxxx yx234( 4 )12 !3 !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)后

3、向O(x2)五點(diǎn)中心O(x4)1iiiyyyx 1iiiyyyx 112iiiyyyx 12432iiiiyyyyx12432iiiiyyyyx 11228812iiiiiyyyyyx精度為O(X2)的高階中心差分算法11221123(4)211242222464iiiiiiiiiiiiiiiyyyyxyyyyyxyyyyyyx精度為O(X4)的高階中心差分算法2112211223211233(4)3211234881216301612813138812395639126iiiiiiiiiiiiiiiiiiiiiiiiiiyyyyyxyyyyyyxyyyyyyyxyyyyyyyyx7.1.2

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

5、-DX(i)= Y(i+2)-2Y(i+1)+Y(i) ， i=1,2 n-2。DX=diff(A,n,dim)：計(jì)算矩陣：計(jì)算矩陣A的的n階差分，階差分，dim=1時(shí)時(shí)(缺省狀態(tài)缺省狀態(tài))，按列計(jì)算差分；按列計(jì)算差分；dim=2，按行計(jì)算差分。，按行計(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=diff(

6、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ù)值梯度，步

7、長為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=gradient(A)Fx = 0 0 0 0 1.00

8、00 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)，對1維向量v以步長h求拉普拉斯算子時(shí)，返回一相同維數(shù)的向

9、量U，且 默認(rèn)的步長為1。 U=4*del2(v,h1,h2)，對矩陣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 數(shù)值積分?jǐn)?shù)值積分7.2.1 數(shù)值積分基本原

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

11、就得到積分精確值。 1( )dlim()nbkkankfxxfx把區(qū)間a，b分割成n等分，像這樣取定步長算積分的方法，稱為定步長積分法法。常見的低階求積分公式復(fù)化矩形公式復(fù)化的梯形公式 復(fù)化的辛普森Simpson公式 10( )bnkkaf x dxS1/21()4()()6kkkkhSf xf xf x()kkSfxh1()()2kkkhSfxfx, (0,1,1)kbaxakh hknn辛普森公式的幾何意義辛普森公式的幾何意義 7.2.2 變步長積分法變步長積分法 計(jì)算積分，可以采取逐步縮小步長計(jì)算積分，可以采取逐步縮小步長h的辦法。的辦法。即先任取步長即先任取步長h進(jìn)行計(jì)算，然后取較小步

12、長進(jìn)行計(jì)算，然后取較小步長 h 進(jìn)進(jìn)行計(jì)算，如果兩次計(jì)算結(jié)果相差較大，則取更小行計(jì)算，如果兩次計(jì)算結(jié)果相差較大，則取更小步長進(jìn)行計(jì)算，如此下去，直到相鄰兩次計(jì)算結(jié)步長進(jìn)行計(jì)算，如此下去，直到相鄰兩次計(jì)算結(jié)果相差不大為止，取最小步長算出的結(jié)果作為積果相差不大為止，取最小步長算出的結(jié)果作為積分值。這種方法稱為變步長積分法。分值。這種方法稱為變步長積分法。 利用兩種步長計(jì)算積分時(shí)，通常取利用兩種步長計(jì)算積分時(shí)，通常取h=h/2 。而每次改變步長后，只需計(jì)算新增節(jié)點(diǎn)處的函數(shù)而每次改變步長后，只需計(jì)算新增節(jié)點(diǎn)處的函數(shù)值，將它們的和乘新步長值，將它們的和乘新步長 。 MATLAB常用的數(shù)值積分函數(shù)參數(shù)名參

13、數(shù)名功能描述功能描述quad一元函數(shù)的數(shù)值積分，采用變步長辛普森Simpson法(低階)quadl一元函數(shù)的數(shù)值積分，采用變步長洛巴托Lobatto法(高階)qualv一元函數(shù)的矢量數(shù)值積分dbquad二重積分（默認(rèn)值為Simpson法）triplequad三重積分（默認(rèn)值為Simpson法） 7.2.3 一元函數(shù)的數(shù)值積分舉例求定積分1. 建立匿名函數(shù)建立匿名函數(shù) f=(x)x./(x.2+1);2. 用用sum函數(shù)實(shí)現(xiàn)復(fù)化矩形法求積函數(shù)實(shí)現(xiàn)復(fù)化矩形法求積1122001log (1)0 .3 4 6 612exd xxxx=0:0.001:1; %步長提高到0.001 y=f(x); sum

14、(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(fun,a,b,tol)其中fun為積分函數(shù)，a,b為積分區(qū)間，tol為積分的誤差閾值，默認(rèn)值為1e-6q

15、uad(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.6827 0.3829 0.2611 0.1974 0.1585矢量積分的結(jié)果是一個(gè)向量，每一個(gè)元素為一個(gè)

16、一元函數(shù)定積分的值。12211exp(),1,2,3,4,522xdx nn7.2.5 二重定積分的數(shù)值求解二重定積分的數(shù)值求解 使用使用MATLAB提供的提供的dblquad函數(shù)就可以直接函數(shù)就可以直接求出上述二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式求出上述二重定積分的數(shù)值解。該函數(shù)的調(diào)用格式為：為： dblquad(f,a,b,c,d,tol,method)該函數(shù)求該函數(shù)求f(x,y)在在a,bc,d區(qū)域上的二重定積分。區(qū)域上的二重定積分。參數(shù)參數(shù)tol，可以采用，可以采用method=quadl的方法使函數(shù)用高的方法使函數(shù)用高階的洛巴托法計(jì)算定積分。階的洛巴托法計(jì)算定積分。例如計(jì)算 f=(x

17、,y)exp(-x.2-y.2)/pi; %歸一化高斯函數(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ù)值解法就是求解所謂的數(shù)值解法就是求解y(x)在在 區(qū)間的近似值區(qū)間的近似值yn的方法，的方法，yn(n=1,2,N)稱為稱為常微分方程的數(shù)值解。自變量常微分方程的數(shù)值解。自變量x的步長一般為定的步長一般為定值值h。0( ,

18、)( )dyfx yaxbdxy ay01Na xxxb7.3.2 常見的數(shù)值方法常見的數(shù)值方法向前歐拉公式向前歐拉公式向后歐拉公式向后歐拉公式梯形公式梯形公式改進(jìn)的歐拉公式改進(jìn)的歐拉公式10(,)0,1( )nnnnyyhfxynNyy a1110(,)0,1( )nnnnyyhfxynNyy a1110(,)(,)0,12( )nnnnnnhyyf xyf xynNyy a1(,)(,)1()2pnnnqnnpnpqyyhfxyyyhfxh yyyy 7.3.3 龍格庫塔法簡介龍格庫塔法簡介 基本思想就是利用在某些點(diǎn)處的值的線性組合基本思想就是利用在某些點(diǎn)處的值的線性組合構(gòu)造公式，使其按泰

19、勒展開后與初值問題的解的泰構(gòu)造公式，使其按泰勒展開后與初值問題的解的泰勒展開式比較，有盡可能多的相同項(xiàng)，從而保證算勒展開式比較，有盡可能多的相同項(xiàng)，從而保證算式有較高的精度。式有較高的精度。常用的四階經(jīng)典的龍格常用的四階經(jīng)典的龍格-庫塔公式庫塔公式112341213243(22)6(,)(,)22(,)22(,)22nnnnnnnnnnhyyKKKKKf x yhhKf xyKhhKf xyKhhKf xyK7.3.4 龍格庫塔法的實(shí)現(xiàn)龍格庫塔法的實(shí)現(xiàn)求解器求解器ODE類類型型特點(diǎn)特點(diǎn)精度精度說明說明ode45非剛性一步算法（只需前一步的結(jié)果），4,5階 Runge-Kutta方法。中大部分場

20、合的首選算法ode23非剛性一步算法，2,3階Runge-Kutta方法。低使用于精度較低的情形ode113非剛性多步法（需要前幾布的結(jié)果），Adams算法。低高計(jì)算時(shí)間比ode45短ode23t適度剛性 采用梯形算法適度剛性情形ode15s剛性多步法，Gears反向數(shù)值積分。低中若ode45失效時(shí)，可嘗試使用ode23s剛性一步法，2階Rosebrock算法精度。低當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短 對于一個(gè)常微分方程組，如果其解相差十分懸殊，就稱之為剛性方程組。對于剛性方程組，為了保持解法的穩(wěn)定，步長選取十分困難，有些解法不再適用。ode函數(shù)調(diào)用格式：t,y = ode ij (ode

21、fun,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.