《數(shù)學(xué)實(shí)驗(yàn) 第4版》課件 4.3 MATLAB數(shù)值積分與微分

第四章數(shù)值分析

實(shí)驗(yàn)4.1插值

實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合數(shù)學(xué)實(shí)驗(yàn)實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分實(shí)驗(yàn)4.4常微分方程的數(shù)值解2實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分一、數(shù)值積分二、數(shù)值微分實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分3一、數(shù)值積分實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分實(shí)驗(yàn)?zāi)康耐ㄟ^本實(shí)驗(yàn)了解數(shù)值積分和數(shù)值微分的方法，會用MATLAB進(jìn)行數(shù)值積分和數(shù)值微分.數(shù)值積分也稱為數(shù)值求積，是求定積分的近似值的數(shù)值方法.數(shù)值積分算法的發(fā)展與完善可以追溯到17世紀(jì)，從最早的牛頓-柯特斯公式到如今的自適應(yīng)積分法和高維積分方法，經(jīng)歷了多個階段不斷的改進(jìn).4實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分

定積分是微積分中的基本計(jì)算方法，但在很多實(shí)際問題中,經(jīng)常會遇到被積函數(shù)的原函數(shù)不能用初等函數(shù)表示；或雖然能找到原函數(shù)但因其很復(fù)雜而難以給出最后的積分結(jié)果；或被積函數(shù)以數(shù)表的形式給出，因此求定積分的數(shù)值解在實(shí)際中應(yīng)用顯得特別重要．?dāng)?shù)學(xué)家們通過不斷的努力和創(chuàng)新，經(jīng)過了幾個世紀(jì)的發(fā)展與完善，提高了數(shù)值積分算法的精度和效率，為科學(xué)計(jì)算和工程應(yīng)用提供了重要的支持.5用數(shù)值方法近似求定積分

的基本思路，就是通過將積分區(qū)間[a,b]劃分成若干小區(qū)間，在每個小區(qū)間上用簡便易求的函數(shù)近似替代被積函數(shù)f(x)，并計(jì)算每個小區(qū)間上的近似函數(shù)所圍成的面積之和來逼近定積分的值.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分如自適應(yīng)辛普森(Simpson)法、自適應(yīng)洛巴托（Lobatto）法、高斯－勒讓德(Gauss-Legendre)法、全局自適應(yīng)求積法等都是經(jīng)常采用的求數(shù)值積分的方法．6MATLAB提供了基于這些算法的相應(yīng)函數(shù)：quad、quadl、quadgk和integral等函數(shù)integral和函數(shù)quad、quadl、quadgk功能基本相同，但前者更強(qiáng)大、更智能化，主要體現(xiàn)在：實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分（1）速度更快；（2）支持積分限為無窮大的積分計(jì)算以及含奇點(diǎn)的積分計(jì)算（quadgk函數(shù)也有此功能）；（3）如果是重積分，integral2和integral3還支持非矩形區(qū)域和非長方體區(qū)域上的積分.71．基于自適應(yīng)求積法的MATLAB實(shí)現(xiàn)基于自適應(yīng)求積法，MATLAB給出了integral函數(shù)來求定積分.該函數(shù)的調(diào)用格式為：I=integral(fun,a,b,Name,Value)fun是函數(shù)句柄．其中a和b分別是定積分的下限和上限．tol用來控制積分精度，默認(rèn)值為tol=0.001．Name,Value是用于指定積分選項(xiàng)的名稱-值對參數(shù)，例如‘AbsTol’和‘RelTol’用于控制絕對和相對誤差容限，默認(rèn)值分別為和.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分I=integral(fun,a,b)或8例10

求定積分解>>I=integral(f10,0,3*pi)↙I=0.9008實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分注

integral函數(shù)也支持無窮區(qū)間，并且能夠處理端點(diǎn)包含奇點(diǎn)的情況.9實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例11

求定積分解>>f11=@(x)exp(-x.^2).*(log(x).^2);I=integral(f11,0,inf)↙I=1.9475注本題積分區(qū)間端點(diǎn)0為奇點(diǎn).102．梯形積分法的MATLAB實(shí)現(xiàn)

在MATLAB中，對于被積函數(shù)以數(shù)表的形式給出的定積分問題用trapz函數(shù)，調(diào)用格式為：其中向量X,Y為等長的兩組向量，定義函數(shù)關(guān)系Y=f(X)．實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分I=trapz(X,Y)一般地，積分區(qū)間是11實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例12已知某次物理實(shí)驗(yàn)測得如下表所示的兩組樣本點(diǎn)：x1.381.562.213.975.517.799.1911.1213.39y3.353.965.128.9811.4617.6324.4129.8332.21現(xiàn)已知變量x和變量y滿足一定的函數(shù)關(guān)系，但此關(guān)系未知，設(shè)y=f(x)，求積分的數(shù)值.解>>X=[1.38,1.56,2.21,3.97,5.51,7.79,9.19,11.12,13.39];Y=[3.35,3.96,5.12,8.98,11.46,17.63,24.41,29.83,32.21];I=trapz(X,Y)↙I=217.103312實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分注函數(shù)關(guān)系式已知的函數(shù)也可以用此命令求定積分的值，需要先生成X，Y的函數(shù)關(guān)系數(shù)據(jù)向量，這種函數(shù)求得的數(shù)值解比函數(shù)integral求得的數(shù)值精確度低.例13

用trapz函數(shù)計(jì)算定積分解>>X=1:0.01:2.5;Y=exp(-X.^2);%生成函數(shù)關(guān)系數(shù)據(jù)向量trapz(X,Y)↙ans=0.1390133.多重積分?jǐn)?shù)值求解的MATLAB實(shí)現(xiàn)使用MATLAB提供的integral2函數(shù)和integral3函數(shù)可以求出矩形區(qū)域上二重積分和長方體區(qū)域上三重積分的的數(shù)值解．這兩個函數(shù)的調(diào)用格式分別為：I=integral2(fun,xmin,xmax,ymin,ymax)，實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)，或I=integral2(fun,xmin,xmax,ymin,ymax,Name,Value)，q=integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax,Name,Value)，14實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分Zmin和zmax分別是z變量的積分下限和上限，這些可以是常數(shù)，也可以是函數(shù)句柄（即可以是x、y的函數(shù)）；ymin和ymax分別是y變量的積分下限和上限，這些可以是常數(shù)，也可以是函數(shù)句柄（即可以是x的函數(shù)）；fun是函數(shù)句柄；其中xmin和xmax分別是x變量的積分下限和上限，必須是有限或無限的實(shí)標(biāo)量值；Name,Value的用法與函數(shù)integral完全相同.15integral2函數(shù)和integral3函數(shù)不僅可以求出矩形區(qū)域上二重積分和長方體區(qū)域上三重積分的數(shù)值解，也可以求出非矩形區(qū)域上二重積分和非長方體區(qū)域上三重積分的數(shù)值解，還支持含奇點(diǎn)的重積分，當(dāng)奇異性位于積分邊界上時，integral2和integral3的性能最佳.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分16實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例14計(jì)算二次積分解>>f14=@(x,y)exp(-x.^2/2).*sin(x.^2+y);I=integral2(f14,-2,2,-1,1)↙I=1.5745注這是函數(shù)

在矩形區(qū)域[-2,2]×[-1,1]上的二重積分.17實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例15計(jì)算二次積分解>>f15=@(x,y)1./(sqrt(x+y).*(1+x+y).^2);ymax=@(x)1-x;%定義y的上限為1-xI=integral2(f15,0,1,0,ymax)↙I=0.2854注這是函數(shù)

在三角形區(qū)域

上的二重積分.積分邊界上含奇點(diǎn)(0,0).18例16計(jì)算三次積分解>>f16=@(x,y,z)y.*sin(x)+z.*cos(x);Q=integral3(f16,0,pi,0,1,-1,1)↙Q=

2.0000注這是函數(shù)在長方體區(qū)域[0,π]×[0,1]×[-1,1]上的三重積分.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分19實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例17計(jì)算三次積分解>>f17=@(x,y,z)1./(1+x+y+z);ymax=@(x)(1-x);%定義y的上限為1-xzmax=@(x,y)(1-x-y);%定義z的上限為1-x-yq=integral3(f17,0,1,0,ymax,0,zmax)↙q=0.0966注這是函數(shù)在四面體區(qū)域

上的三重積分.積分邊界上含奇點(diǎn)(0,0,0).20二、數(shù)值微分實(shí)際中常遇到僅給出了一系列離散點(diǎn)及相應(yīng)函數(shù)值的列表型函數(shù)的求導(dǎo)問題，這就需要用這些離散點(diǎn)的函數(shù)值推算函數(shù)在某點(diǎn)的導(dǎo)數(shù)或高階導(dǎo)數(shù)的近似值，這種方法稱為數(shù)值微分.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分對于難以求導(dǎo)的復(fù)雜函數(shù)，也可以用數(shù)值微分求導(dǎo)，不過需要先由函數(shù)表達(dá)式生成離散的數(shù)據(jù)列表.通常用以下三種思路建立數(shù)值微分公式：21實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分(1)差商近似數(shù)值微分：從導(dǎo)數(shù)定義出發(fā)，通過近似處理，得到數(shù)值微分；(2)插值型數(shù)值微分：利用本章實(shí)驗(yàn)4.1介紹的插值公式得到近似代替該函數(shù)的較簡單函數(shù)，對其求導(dǎo)得到要求導(dǎo)數(shù)的近似值.(3)擬合型數(shù)值微分：利用本章實(shí)驗(yàn)4.2介紹的數(shù)據(jù)擬合的方法得到近似代替該函數(shù)的較簡單函數(shù)，對其求導(dǎo)得到要求導(dǎo)數(shù)的近似值.22實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分1．差商近似數(shù)值微分（1）數(shù)值微分與微商導(dǎo)數(shù)定義為假設(shè)h>0，引進(jìn)記號函數(shù)f

(x)在x點(diǎn)處以h

為步長的向前差分函數(shù)f

(x)在x點(diǎn)處以h

為步長的向前差商當(dāng)步長h足夠小時，有稱為向前差商數(shù)值微分公式23實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分類似可得向后差商數(shù)值微分公式中心差商數(shù)值微分公式其中中心差商公式精度較高.24DX=diff(X)：計(jì)算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1DX=diff(X,n)：計(jì)算X的n階向前差分．例如，diff(X,2)=diff(diff(X))DX=diff(A,n,dim)：計(jì)算矩陣A的n階差分，dim=1時(默認(rèn)狀態(tài))，按列

計(jì)算差分；dim=2，按行計(jì)算差分．實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分（2）差分的MATLAB實(shí)現(xiàn)在MATLAB中，沒有直接提供求數(shù)值導(dǎo)數(shù)的函數(shù)，只有計(jì)算向前差分的函數(shù)diff，利用它可以求出微商，從而得到要求的近似導(dǎo)數(shù).diff的調(diào)用格式為:25例18根據(jù)下表所示年份出生人口數(shù)，計(jì)算出生人口年增長率解差商近似求數(shù)值微分的程序如下：實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分年份193019351940194519501955196019651970人口/萬650781914100514711861146824792801年份197519801985199019952000200520102015人口/萬211418392043262116931379161715741655>>t=[1930:5:2015];p=[650781914100514711861146824792801211418392043262116931379161715741655];26dt=diff(t);%求時間t的差分dp=diff(p);%求人口p的差分q=dp./dt↙%利用差商求數(shù)值導(dǎo)數(shù),即出生人口增長率列1至626.200026.600018.200093.200078.0000-78.6000列7至12202.200064.4000-137.4000-55.000040.8000115.6000列13至17-185.6000-62.800047.6000-8.600016.2000差商近似法是最簡單的數(shù)值微分方法，在實(shí)際中十分常用，但其精確度不高,誤差較大.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分27實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分2．插值型數(shù)值微分插值型數(shù)值微分是差商近似法的推廣.當(dāng)函數(shù)可微性不太好時，利用樣條插值進(jìn)行數(shù)值微分要比多項(xiàng)式插值更適宜.僅就三次樣條插值方法說明數(shù)值微分過程：離散數(shù)據(jù)三次樣條插值函數(shù)pp的導(dǎo)數(shù)pppp在點(diǎn)xi的導(dǎo)數(shù)值fnder是對樣條函數(shù)求導(dǎo)，fnval用來計(jì)算樣條函數(shù)的函數(shù)值.28實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例19某液體冷卻時，溫度隨時間的變化數(shù)據(jù)如下表所示.試分別計(jì)算t=2,3,4min及t=1.5,2.5,4.5min時的降溫速率.t/min012345T/oC92.085.379.574.570.267.0分析前者是計(jì)算節(jié)點(diǎn)處的一階導(dǎo)數(shù)，后者是計(jì)算非節(jié)點(diǎn)處

的一階導(dǎo)數(shù).解三次樣條插值函數(shù)求數(shù)值微分的程序如下：>>t=[0:5];T=[92.0,85.3,79.5,74.5,70.2,67.0];29實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分p=spline(t,T);%生成三次樣條插值函數(shù)pp=fnder(p);%生成三次樣條插值函數(shù)的導(dǎo)函數(shù)t1=[2,3,4,1.5,2.5,4.5];dT=fnval(pp,t1);%計(jì)算導(dǎo)函數(shù)在t1處的導(dǎo)數(shù)值disp('相應(yīng)時間時的降溫速率：')disp([t1;dT])↙相應(yīng)時間時的降溫速率：2.00003.00004.00001.50002.50004.5000-5.3722-4.6722-3.8389-5.7972-4.9889-3.2222注插值型數(shù)值微分不但適用于求節(jié)點(diǎn)處的導(dǎo)數(shù)，還可以求非節(jié)點(diǎn)處的導(dǎo)數(shù).30實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分3．?dāng)M合型數(shù)值微分如果離散點(diǎn)上的數(shù)據(jù)有不容忽視的隨機(jī)誤差，應(yīng)該用曲線擬合代替函數(shù)插值，然后用擬合曲線的導(dǎo)數(shù)作為所求導(dǎo)數(shù)的近似值，這種做法可以起到減少隨機(jī)誤差的作用.僅就多項(xiàng)式擬合方法說明數(shù)值微分過程:離散數(shù)據(jù)多項(xiàng)式擬合函數(shù)導(dǎo)函數(shù)pppp在點(diǎn)xi的導(dǎo)數(shù)值polyder是對多項(xiàng)式求導(dǎo).31實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例20一底面面積為常數(shù)

S的正圓柱體水塔的某一天

0~9

點(diǎn)的水位測量記錄(這一時段沒有給水塔充水)如下表所示,根據(jù)該表估計(jì)這一時段任何時刻從水塔流出的水流量.分析由于水塔截面積是常數(shù)，為簡單起見，計(jì)算中將流量定義為單位時間流出水的高度，即水位對時間變化率的絕對值(水位是下降的).時間/h00.921.842.953.874.985.907.017.938.97水位/cm96894893191389888186985283982232實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分解方法一多項(xiàng)式擬合求數(shù)值微分的程序如下：>>t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97];h=[968,948,931,913,898,881,869,852,839,822];A=polyfit(t,h,3);%3次多項(xiàng)式擬合B=polyder(A);%對擬合多項(xiàng)式求導(dǎo)tp=0:0.1:9;x=-polyval(B,tp)↙%求tp時刻的水流量x=列1至722.107921.838521.573921.313921.058720.808220.5623列8至1420.321220.084919.853219.626219.404019.186418.973633實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分列

15至2118.765518.562118.363418.169517.980217.795717.6158列22至2817.440717.270317.104616.943616.787316.635816.4889列29至3516.346816.209416.076715.948715.825415.706815.5929列36至4215.483815.379415.279615.184615.094315.008814.9279列43至4914.851714.780314.713514.651514.594214.541614.4937列50至5614.450614.412114.378414.349314.325014.305414.290534實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分列

57至6314.280314.274814.274114.278014.286714.300114.3182列64至7014.341014.368514.400714.437714.479314.525714.5768列71至7714.632614.693114.758314.828214.902914.982215.0663列78至8415.155115.248515.346715.449715.557315.669615.7867列85至9115.908516.034916.166116.302016.442716.588016.7380我們可以用給定時段的用水量

968-822=146檢驗(yàn)計(jì)算結(jié)果：35實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分在上述程序下繼續(xù)運(yùn)行>>y=0.1*trapz(x)%用數(shù)值積分計(jì)算給定時段的總用水量，

積分步長為0.1↙y=146.1815與用水量的絕對誤差為146-146.1815=0.1815.36實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分方法二差商近似求數(shù)值微分的程序如下：>>t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97];h=[968,948,931,913,898,881,869,852,839,822];dt=diff(t);%求時間t的差分dh=diff(h);%求水位h的差分q=dh./dt↙%利用差分求數(shù)值導(dǎo)數(shù)，即得已知時刻水流量q=-21.7391-18.4783-16.2162-16.3043-15.3153-13.0435-15.3153-14.1304-16.3462>>u=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93];A=polyfit(u,q,3);%用3次多項(xiàng)式擬合流量函數(shù)的系數(shù)t=0:0.1:9;37實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分s=-polyval(A,t)↙%輸出t時刻的水流量值s=列1至721.339821.048520.763920.486020.214819.950219.6921列8至1419.440619.195618.957118.724918.499218.279818.0666列15至2117.859717.659117.464617.276217.094016.917816.7476列22至2816.583316.425116.272716.126115.985415.850415.7212列29至3515.597615.479715.367515.260815.159615.063914.973738實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分列36至4214.888814.809414.735314.666414.602814.544414.4912列43至4914.443114.400214.362214.329314.301314.278314.2601列50至5614.246814.238314.234614.235514.241214.251514.2665列57至6314.286014.310014.338514.371414.408814.4