1第四章數(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.1插值一、拉格朗日（Lagrange）插值法二、分段線性插值三、三次樣條插值四、應(yīng)用舉例3實(shí)驗(yàn)4.1插值一、拉格朗日（Lagrange）插值法實(shí)際應(yīng)用中，經(jīng)常遇到下面的問題：通過實(shí)驗(yàn)、測(cè)量等方法得到一些離散數(shù)據(jù)，這些數(shù)據(jù)從數(shù)學(xué)的角度可以看作是由某個(gè)函數(shù)產(chǎn)生的，如何利用這些數(shù)據(jù)來尋找這個(gè)函數(shù)滿足要求精度，且相對(duì)簡(jiǎn)單的近似表達(dá)式呢？我們介紹以下三種方法：多項(xiàng)式插值、分段線性插值、三次樣條插值.4這時(shí)我們稱為插值函數(shù),為插值節(jié)點(diǎn).作為的近似表達(dá)式,使?jié)M足的一組測(cè)量數(shù)據(jù)設(shè)函數(shù),要尋求一個(gè)函數(shù)1.拉格朗日插值多項(xiàng)式在我們所學(xué)的函數(shù)類型中，多項(xiàng)式相對(duì)比較簡(jiǎn)單，用多項(xiàng)式作為插值函數(shù)是常用的方法，也稱為多項(xiàng)式插值法.實(shí)驗(yàn)4.1插值5利用上述方法求n次插值多項(xiàng)式計(jì)算比較麻煩，由插值多項(xiàng)式的唯一性，無論用什么方法，由n+1個(gè)節(jié)點(diǎn)確定的n次插值多項(xiàng)式都是相同的.當(dāng)n+1個(gè)節(jié)點(diǎn)互不相同時(shí)，由線性代數(shù)的知識(shí)，方程組(3)唯一確定

的一組系數(shù)

。由此可見,n+1個(gè)節(jié)點(diǎn)可以唯一確定一個(gè)n次插值多項(xiàng)式.實(shí)驗(yàn)4.1插值6比較方便的方法是構(gòu)造一組基函數(shù)：（4）(5)令（6）就是滿足條件(1)的n次插值多項(xiàng)式,稱為n次Lagrange(拉格朗表示，即日)插值多項(xiàng)式，通常用（7）實(shí)驗(yàn)4.1插值72.多項(xiàng)式插值的誤差估計(jì)關(guān)于可以證明有下列結(jié)論成立.(8)若

在

上連續(xù)，

在

內(nèi)存在，若

是

上的

個(gè)互異節(jié)點(diǎn)，則插值多項(xiàng)式

，對(duì)任意點(diǎn)

，插值余項(xiàng)為其中

，且與

有關(guān)定理(9)實(shí)驗(yàn)4.1插值8實(shí)驗(yàn)4.1插值functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endy3.拉格朗日插值法的MATLAB實(shí)現(xiàn)把這個(gè)程序存盤起來，作拉格朗日插值計(jì)算時(shí)可以隨時(shí)調(diào)用.編寫拉格朗日插值法程序：注意:程序中n表示節(jié)點(diǎn)的個(gè)數(shù)，以數(shù)組x0,y0輸入，m表示插值點(diǎn)的個(gè)數(shù)，以數(shù)組x輸入.輸出數(shù)組y為m個(gè)插值.9分別用n=8、10的等距分點(diǎn)進(jìn)行多項(xiàng)式插值，例1繪制f(x)及插值多項(xiàng)式的圖形.解x=[-1:0.25:1];y=1./(1+9*x.^2);x0=[-1:0.05:1];y0=lagrange(x,y,x0);plot(x0,y0,'r--')holdon,y1=1./(1+9*x0.^2);plot(x0,y1),holdoff↙在命令窗口輸入：取區(qū)間[-1,1]的8等分點(diǎn)作為節(jié)點(diǎn),實(shí)驗(yàn)4.1插值函數(shù)f(x)及插值多項(xiàng)式的圖形如圖4.1所示.10實(shí)驗(yàn)4.1插值圖4.111取區(qū)間[-1,1]的10等分點(diǎn)作為節(jié)點(diǎn),類似地可得f(x)及插值多項(xiàng)式,如圖4.2所示.實(shí)驗(yàn)4.1插值圖4.212一般總認(rèn)為插值多項(xiàng)式的次數(shù)越高,逼近的精度越高.但事實(shí)上并非如此,由圖4.1和圖4.2可以看出，離零越遠(yuǎn),近似效果越差,以致完全失真.這一現(xiàn)象被稱為Runge現(xiàn)象.Runge現(xiàn)象表明,盲目采用提高插值多項(xiàng)式的次數(shù)的方法來減少誤差是不可取的.因此,我們只有用縮短插值區(qū)間,分段進(jìn)行插值來達(dá)到減少誤差的目的.與f（x）的近似程度比較好，但當(dāng)x實(shí)驗(yàn)4.1插值13實(shí)驗(yàn)4.1插值二、分段線性插值

分段線性插值,就是用連接彼此相鄰兩節(jié)點(diǎn)的直線段形成的折線作為插值函數(shù).MATLAB對(duì)分段線性插值提供了插值函數(shù),見下表.函數(shù)功能yi=interp1(x,y,xi)對(duì)節(jié)點(diǎn)(x,y)插值，求xi處的內(nèi)插值.yi=interp1(x,y,xi,'method')采用指定的方法，對(duì)節(jié)點(diǎn)(x,y)插值，求xi處的內(nèi)插值.14實(shí)驗(yàn)4.1插值用n=20的等距分點(diǎn)進(jìn)行線插值，繪制f(x)及插值例2多項(xiàng)式的圖形.在命令窗口輸入：x=-1:0.1:1;y=1./(1+9*x.^2);xi=-1:0.1:1;yi=interp1(x,y,xi);plot(x,y,'r-',xi,yi,'*')↙解利用plot函數(shù)檢驗(yàn)插值效果，由圖4.3可見，在[-1,1]內(nèi)線性插值函數(shù)收斂于函數(shù)

f(x).15

一般地，線性插值函數(shù)都有良好的收斂性.數(shù)學(xué)、物理中用的特殊函數(shù)表,計(jì)算機(jī)繪圖都采用了分段線性插值的原理.實(shí)驗(yàn)4.1插值圖4.316三、三次樣條插值

分段線性插值雖然有良好的收斂性,但是由于在節(jié)點(diǎn)處不光滑,實(shí)用上受到了一定的限制.下面介紹一種在實(shí)際應(yīng)用中,常用的提高光滑度的方法:三次樣條插值.(1)在每一個(gè)小區(qū)間上是三次多項(xiàng)式且則稱S(x)為三次樣條插值函數(shù).若函數(shù)S(x)滿足下列條件：給定[a,b]上的n+1個(gè)節(jié)點(diǎn)滿足實(shí)驗(yàn)4.1插值17實(shí)驗(yàn)4.1插值其中為待定常數(shù).由條件（2）共4n-2個(gè)條件，還需要再給兩個(gè)條件，才能唯一確定一個(gè)三次樣條插值函數(shù).最常用的是增設(shè)邊界條件:也稱為自然條件.18實(shí)驗(yàn)4.1插值在MATLAB中三次樣條插值函數(shù)為：yi=interp1(x,y,xi,'spline'),

或yi=spline(x,y,xi)其中變量的含義與分段線性插值相同.19例3在一天24小時(shí)內(nèi),從零點(diǎn)開始每間隔2小時(shí)測(cè)得的環(huán)境溫度為(攝氏度)12,9,9,10,18,24,28,27,25,20,18,15,13推測(cè)在每1s時(shí)的溫度，并描繪溫度曲線.在命令窗口輸入：解t=0:2:24;T=[129910182428272520181513];plot(t,T,'*')ti=0:1/3600:24;T1i=interp1(t,T,ti);holdon,plot(ti,T1i)T2i=interp1(t,T,ti,'spline');holdon,plot(ti,T2i,'r-')↙實(shí)驗(yàn)4.1插值20輸入不同的時(shí)刻ti便可得到相應(yīng)的溫度值Ti.利用線性插值方法描繪的溫度曲線，如圖4.4所示.通過圖上兩條曲線比較可以看出,三次樣條插值函數(shù)在整個(gè)區(qū)間上有較好的收斂性;光滑性比分段線性插值有較大的提高,因此應(yīng)用比較廣泛.缺點(diǎn)是誤差估計(jì)比較困難.這正像我們的人生在發(fā)揚(yáng)優(yōu)點(diǎn)修正錯(cuò)誤中不斷砥礪前行.圖4.4實(shí)驗(yàn)4.1插值21四、應(yīng)用舉例數(shù)據(jù)加細(xì)問題在機(jī)械制造加工中,經(jīng)常遇到數(shù)據(jù)加細(xì)的問題.例如，在現(xiàn)代機(jī)械工業(yè)中，用計(jì)算機(jī)程序控制加工機(jī)器零件時(shí)，根據(jù)設(shè)計(jì)可以給出零件外形曲線上某些點(diǎn)，加工時(shí)為控制每步刀的走向及步長(zhǎng)就要算出零件外形曲線上其他點(diǎn)的函數(shù)值，加工出外表光滑的零件，這就涉及到數(shù)據(jù)加細(xì)的問題.例4

在飛機(jī)的機(jī)翼加工時(shí)，由于機(jī)翼尺寸很大，通常在圖紙上只能標(biāo)出部分關(guān)鍵點(diǎn)的數(shù)據(jù).某型號(hào)飛機(jī)的機(jī)翼上緣輪廓線的部分?jǐn)?shù)據(jù)如下：實(shí)驗(yàn)4.1插值x0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00y0.005.238.1011.9716.1517.1016.3414.6312.166.697.033.990.00對(duì)表中的數(shù)據(jù)進(jìn)行細(xì)化，并畫出機(jī)翼的上輪廓線.22解輸入：x=[0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00];y=[0.005.238.1011.9716.1517.1016.3414.6312.169.697.033.990.00];xi=[0:0.001:190];yi=interp1(x,y,xi,'spline');plot(xi,yi)↙實(shí)驗(yàn)4.1插值機(jī)翼的上輪廓線，如圖所示.圖4.523實(shí)驗(yàn)4.1插值例5

天文學(xué)家在1914年8月份的7次觀測(cè)中，測(cè)得地球與金星之間距離（單位：米），并取其常用對(duì)數(shù)值與日期的一組歷史數(shù)據(jù)如下所示,試推斷何時(shí)金星與地球的距離(米)的對(duì)數(shù)值為9.9352.日期:18202224262830距離對(duì)數(shù):9.96189.95449.94689.93919.93129.92329.9150解由于對(duì)數(shù)值9.9352位于24和26兩天所對(duì)應(yīng)的對(duì)數(shù)值之間,所以對(duì)上述數(shù)據(jù)用三次樣條插值加細(xì)為步長(zhǎng)為1的數(shù)據(jù),輸入:x=[18:2:30];y=[9.96189.95449.94689.93919.93129.92329.9150];xi=[18:1:30];yi=interp1(x,y,xi,'spline')；A=[xi;yi]↙24經(jīng)三次樣條插值推斷，25日時(shí)金星與地球的距離(米)的對(duì)數(shù)值為9.9352.A=18.000019.000020.000021.000022.000023.000024.00009.96189.95819.95449.95069.94689.94309.939125.000026.000027.000028.000029.000030.00009.93529.93129.92729.92329.91919.9150實(shí)驗(yàn)4.1插值25第四章數(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ù)值解26實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合一、離散數(shù)據(jù)的多項(xiàng)式擬合二、曲線擬合的線性最小二乘法三、應(yīng)用舉例數(shù)學(xué)實(shí)驗(yàn)27實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合一、離散數(shù)據(jù)的多項(xiàng)式擬合p=polyfit(x,y,n)用多項(xiàng)式擬合一組離散數(shù)據(jù)就是尋找一組多項(xiàng)式的系數(shù)使得多項(xiàng)式能夠較好的擬合這組數(shù)據(jù).它與實(shí)驗(yàn)4.1的插值法不同,數(shù)據(jù)不能保證都在擬合多項(xiàng)式曲線上，但能使整體擬合誤差較小.在MATLAB中,多項(xiàng)式擬合可以通過polyfit函數(shù)來實(shí)現(xiàn),該函數(shù)的調(diào)用格式為p是多項(xiàng)式系數(shù)按降冪排列得出的行向量.28利用poly2sym(p)可以得出相應(yīng)多項(xiàng)式的表達(dá)式.如果要計(jì)算擬合多項(xiàng)式在x處的值y,輸入y=polyval(p,x),就可輸出y的值.求本章例4中，數(shù)據(jù)組的3次、6次和8次多項(xiàng)式并作圖例6解輸入：x=[0.004.749.0519.0038.0057.0076.0095.00114.00133.00152.00171.00190.00];y=[0.005.238.1011.9716.1517.1016.3414.6312.166.697.033.990.00];p3=polyfit(x,y,3);y1=polyval(p3,x);p6=polyfit(x,y,6);y2=polyval(p6,x);plot(x,y,'*',x,y1,'--',x,y2,'-.')↙實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合29從圖4.6可見，6次擬合多項(xiàng)式比3次多項(xiàng)式的擬合效果好.實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合30在命令窗口輸入：P6=polyfit(x,y,6);y2=polyval(p6,x);P8=polyfit(x,y,8);y3=polyval(p8,x);plot(x,y,'*'x,y2,'—.',x,y3)↙實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合由圖4.7所示，8次多項(xiàng)式比6次的擬合效果更好.可見，隨著多項(xiàng)式次數(shù)的不斷增加，擬合的效果也越來越好，當(dāng)擬合多項(xiàng)式的次數(shù)就能得出較好的效果.那么利用多項(xiàng)式進(jìn)行擬合時(shí)，是否多項(xiàng)式的次數(shù)越高擬合效果一定就越好呢？31設(shè)已知數(shù)據(jù)來自函數(shù)8次多項(xiàng)式擬合，并作圖.例7解試用生成的數(shù)據(jù)進(jìn)行3次、6次和在命令窗口輸入：x0=-1:.01:1;y0=1./(1+9*x0.^2);p3=polyfit(x0,y0,3);y1=polyval(p3,x0);p6=polyfit(x0,y0,6);y2=polyval(p6,x0);p8=polyfit(x0,y0,8);y3=polyval(p8,x0);plot(x0,y0,'*',x0,y1,'--',x0,y2,'-.',x0,y3)↙32由圖4.8可以看出,多項(xiàng)式擬合的效果并不一定總是很精確的.下面我們來介紹另一種方法—曲線擬合的線性最小二乘法.實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合圖4.833二、曲線擬合的線性最小二乘法最小,就是曲線擬合得最好.曲線擬合常用的方法是線性最小二乘法.

曲線擬合（curvefitting）是指選擇適當(dāng)?shù)那€類型來擬合觀測(cè)數(shù)據(jù)，并用擬合的曲線方程分析兩變量間的關(guān)系.實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合已知某函數(shù)的一組測(cè)量數(shù)據(jù)根據(jù)這組數(shù)據(jù)尋求曲線逼近曲線因?yàn)闇y(cè)量時(shí)可能產(chǎn)生誤差,所以我們不要求都經(jīng)過這些點(diǎn),只要與的距離最為接近,即34一般情況下,可以假設(shè)數(shù)據(jù)擬合曲線為將（2）代入(1),上述問題轉(zhuǎn)化為根據(jù)多元函數(shù)極值的必要條件實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合35可以證明方程組(4)的系數(shù)矩陣是可逆的,則方程組(4)有唯一解于是有這種求擬合曲線的方法稱為曲線擬合的線性最小二乘法.實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合36在MATLAB優(yōu)化工具箱中，提供了函數(shù)lsqcurvefit,求解最小二乘曲線擬合問題.該函數(shù)的調(diào)用格式為a=lsqcurvefit(fun,a0,x,y)其中，fun是自定義函數(shù)的MATLAB表示，可以用inline('函數(shù)內(nèi)容'，自變量列表)；a0是a的初始預(yù)測(cè)值，使用方法見例8.實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合37三、應(yīng)用舉例用切削機(jī)床進(jìn)行金屬品加工時(shí)，為了適當(dāng)?shù)恼{(diào)整機(jī)床，需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度，得數(shù)據(jù)如下表所示：例8切削時(shí)間t/h012345678刀具厚度y/cm30.029.1

28.428.1

28.0

27.727.5

27.227.0切削時(shí)間t/h910111213141516刀具厚度y/cm26.8

26.5

26.326.125.725.324.824.0假設(shè)經(jīng)驗(yàn)公式是試用最小二乘法確定實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合38解定義函數(shù)f，在命令窗口輸入:f=inline('a(1).*t.^3+a(2).*t.^2+a(3).*t+a(4)','a','t')確定經(jīng)驗(yàn)公式的系數(shù)并作圖，在命令窗口輸入:t=[0:1:16];y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0];a0=[0,0,0,0]a=lsqcurvefit(f,a0,t,y)y1=a(1).*t.^3+a(2).*t.^2+a(3).*t+a(4);plot(t,y,'*',t,y1)↙實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合a=-0.00290.0678-0.713329.824939系數(shù)散點(diǎn)及擬合曲線圖，實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合如圖4.9所示.40例9在加壓實(shí)驗(yàn)中的應(yīng)力-應(yīng)變關(guān)系測(cè)試點(diǎn)的數(shù)據(jù)如下表所示：已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描述.即假設(shè)使用最小二乘法確定參數(shù)實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合41解選取指數(shù)函數(shù)作擬合時(shí),在擬合前需作變量代換,化為將(6)變形后取對(duì)數(shù)，得令式（7）化為在命令窗口輸入：x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6];y=[3.100*1.0e+32.470*1.0e+31.953*1.0e+31.515*1.0e+31.217*1.0e+3];實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合42w=[1.552.472.933.032.89];plot(x,w,'*')↙holdon,y1=exp(8.3020)*x.*exp(-495.4888*x);plot(x,y1,'r-'),holdoff↙實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合z=log(y)↙z=8.03927.81207.57717.32327.1041a=polyfit(x,z,1)↙a=-495.48888.3020k1=exp(8.2959)↙k1=4.0319e+03所求的擬合曲線(見圖4.10)所以43指數(shù)函數(shù)經(jīng)常應(yīng)用在預(yù)應(yīng)力混凝土梁的分析中，作為應(yīng)力－應(yīng)變關(guān)系的數(shù)學(xué)模型.圖4-10實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合44在實(shí)際應(yīng)用中常見的擬合曲線有直線多項(xiàng)式(一般n=2,3,不宜過高)雙曲線（一支）指數(shù)曲線實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合第四章數(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ù)值解46實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分一、數(shù)值積分二、數(shù)值微分實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分47一、數(shù)值積分實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分實(shí)驗(yàn)?zāi)康耐ㄟ^本實(shí)驗(yàn)了解數(shù)值積分和數(shù)值微分的方法，會(huì)用MATLAB進(jìn)行數(shù)值積分和數(shù)值微分.數(shù)值積分也稱為數(shù)值求積，是求定積分的近似值的數(shù)值方法.數(shù)值積分算法的發(fā)展與完善可以追溯到17世紀(jì)，從最早的牛頓-柯特斯公式到如今的自適應(yīng)積分法和高維積分方法，經(jīng)歷了多個(gè)階段不斷的改進(jìn).48實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分

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

的基本思路，就是通過將積分區(qū)間[a,b]劃分成若干小區(qū)間，在每個(gè)小區(qū)間上用簡(jiǎn)便易求的函數(shù)近似替代被積函數(shù)f(x)，并計(jì)算每個(gè)小區(qū)間上的近似函數(shù)所圍成的面積之和來逼近定積分的值.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分如自適應(yīng)辛普森(Simpson)法、自適應(yīng)洛巴托（Lobatto）法、高斯－勒讓德(Gauss-Legendre)法、全局自適應(yīng)求積法等都是經(jīng)常采用的求數(shù)值積分的方法．50MATLAB提供了基于這些算法的相應(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ū)域和非長(zhǎng)方體區(qū)域上的積分.511．基于自適應(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)的名稱-值對(duì)參數(shù)，例如‘AbsTol’和‘RelTol’用于控制絕對(duì)和相對(duì)誤差容限，默認(rèn)值分別為和.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分I=integral(fun,a,b)或52例10

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

integral函數(shù)也支持無窮區(qū)間，并且能夠處理端點(diǎn)包含奇點(diǎn)的情況.53實(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).542．梯形積分法的MATLAB實(shí)現(xiàn)

在MATLAB中，對(duì)于被積函數(shù)以數(shù)表的形式給出的定積分問題用trapz函數(shù)，調(diào)用格式為：其中向量X,Y為等長(zhǎng)的兩組向量，定義函數(shù)關(guān)系Y=f(X)．實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分I=trapz(X,Y)一般地，積分區(qū)間是55實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例12已知某次物理實(shí)驗(yàn)測(cè)得如下表所示的兩組樣本點(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.103356實(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.1390573.多重積分?jǐn)?shù)值求解的MATLAB實(shí)現(xiàn)使用MATLAB提供的integral2函數(shù)和integral3函數(shù)可以求出矩形區(qū)域上二重積分和長(zhǎng)方體區(qū)域上三重積分的的數(shù)值解．這兩個(gè)函數(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)，58實(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完全相同.59integral2函數(shù)和integral3函數(shù)不僅可以求出矩形區(qū)域上二重積分和長(zhǎng)方體區(qū)域上三重積分的數(shù)值解，也可以求出非矩形區(qū)域上二重積分和非長(zhǎng)方體區(qū)域上三重積分的數(shù)值解，還支持含奇點(diǎn)的重積分，當(dāng)奇異性位于積分邊界上時(shí)，integral2和integral3的性能最佳.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分60實(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]上的二重積分.61實(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).62例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ù)在長(zhǎng)方體區(qū)域[0,π]×[0,1]×[-1,1]上的三重積分.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分63實(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).64二、數(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ù)值積分與微分對(duì)于難以求導(dǎo)的復(fù)雜函數(shù)，也可以用數(shù)值微分求導(dǎo)，不過需要先由函數(shù)表達(dá)式生成離散的數(shù)據(jù)列表.通常用以下三種思路建立數(shù)值微分公式：65實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分(1)差商近似數(shù)值微分：從導(dǎo)數(shù)定義出發(fā)，通過近似處理，得到數(shù)值微分；(2)插值型數(shù)值微分：利用本章實(shí)驗(yàn)4.1介紹的插值公式得到近似代替該函數(shù)的較簡(jiǎn)單函數(shù)，對(duì)其求導(dǎo)得到要求導(dǎo)數(shù)的近似值.(3)擬合型數(shù)值微分：利用本章實(shí)驗(yàn)4.2介紹的數(shù)據(jù)擬合的方法得到近似代替該函數(shù)的較簡(jiǎn)單函數(shù)，對(duì)其求導(dǎo)得到要求導(dǎo)數(shù)的近似值.66實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分1．差商近似數(shù)值微分（1）數(shù)值微分與微商導(dǎo)數(shù)定義為假設(shè)h>0，引進(jìn)記號(hào)函數(shù)f

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

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

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

為步長(zhǎng)的向前差商當(dāng)步長(zhǎng)h足夠小時(shí)，有稱為向前差商數(shù)值微分公式67實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分類似可得向后差商數(shù)值微分公式中心差商數(shù)值微分公式其中中心差商公式精度較高.68DX=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時(shí)(默認(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)用格式為:69例18根據(jù)下表所示年份出生人口數(shù)，計(jì)算出生人口年增長(zhǎng)率解差商近似求數(shù)值微分的程序如下：實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分年份193019351940194519501955196019651970人口/萬650781914100514711861146824792801年份197519801985199019952000200520102015人口/萬211418392043262116931379161715741655>>t=[1930:5:2015];p=[650781914100514711861146824792801211418392043262116931379161715741655];70dt=diff(t);%求時(shí)間t的差分dp=diff(p);%求人口p的差分q=dp./dt↙%利用差商求數(shù)值導(dǎo)數(shù),即出生人口增長(zhǎng)率列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差商近似法是最簡(jiǎn)單的數(shù)值微分方法，在實(shí)際中十分常用，但其精確度不高,誤差較大.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分71實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分2．插值型數(shù)值微分插值型數(shù)值微分是差商近似法的推廣.當(dāng)函數(shù)可微性不太好時(shí)，利用樣條插值進(jìn)行數(shù)值微分要比多項(xiàng)式插值更適宜.僅就三次樣條插值方法說明數(shù)值微分過程：離散數(shù)據(jù)三次樣條插值函數(shù)pp的導(dǎo)數(shù)pppp在點(diǎn)xi的導(dǎo)數(shù)值fnder是對(duì)樣條函數(shù)求導(dǎo)，fnval用來計(jì)算樣條函數(shù)的函數(shù)值.72實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例19某液體冷卻時(shí)，溫度隨時(shí)間的變化數(shù)據(jù)如下表所示.試分別計(jì)算t=2,3,4min及t=1.5,2.5,4.5min時(shí)的降溫速率.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];73實(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)時(shí)間時(shí)的降溫速率：')disp([t1;dT])↙相應(yīng)時(shí)間時(shí)的降溫速率：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ù).74實(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是對(duì)多項(xiàng)式求導(dǎo).75實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分例20一底面面積為常數(shù)S的正圓柱體水塔的某一天

0~9

點(diǎn)的水位測(cè)量記錄(這一時(shí)段沒有給水塔充水)如下表所示,根據(jù)該表估計(jì)這一時(shí)段任何時(shí)刻從水塔流出的水流量.分析由于水塔截面積是常數(shù)，為簡(jiǎn)單起見，計(jì)算中將流量定義為單位時(shí)間流出水的高度，即水位對(duì)時(shí)間變化率的絕對(duì)值(水位是下降的).時(shí)間/h00.921.842.953.874.985.907.017.938.97水位/cm96894893191389888186985283982276實(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);%對(duì)擬合多項(xiàng)式求導(dǎo)tp=0:0.1:9;x=-polyval(B,tp)↙%求tp時(shí)刻的水流量x=列1至722.107921.838521.573921.313921.058720.808220.5623列8至1420.321220.084919.853219.626219.404019.186418.973677實(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.290578實(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我們可以用給定時(shí)段的用水量968-822=146檢驗(yàn)計(jì)算結(jié)果：79實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分在上述程序下繼續(xù)運(yùn)行>>y=0.1*trapz(x)%用數(shù)值積分計(jì)算給定時(shí)段的總用水量，

積分步長(zhǎng)為0.1↙y=146.1815與用水量的絕對(duì)誤差為146-146.1815=0.1815.80實(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);%求時(shí)間t的差分dh=diff(h);%求水位h的差分q=dh./dt↙%利用差分求數(shù)值導(dǎo)數(shù)，即得已知時(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;81實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分s=-polyval(A,t)↙%輸出t時(shí)刻的水流量值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.973782實(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.450514.4966列64至7014.546914.601514.660214.723214.790214.861414.9366列71至7715.015815.098915.186015.277015.371815.470415.572783實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分列78至8415.678815.788615.902016.019016.139616.263716.3913列85至9116.522316.656716.794416.935517.079917.227517.3783我們?nèi)匀豢梢詸z驗(yàn)結(jié)果：w=0.1*trapz(s)%用數(shù)值積分計(jì)算給定時(shí)段的總用水量，

積分步長(zhǎng)為0.1↙w=144.0565與用水量的絕對(duì)誤差為146-144.0565=1.9435.84實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分注比較兩種方法發(fā)現(xiàn)，擬合法誤差更小，更精確.擬合法的不足之處在于：當(dāng)給出等距節(jié)點(diǎn)時(shí)，不如差商近似法通用性強(qiáng)，擬合的多項(xiàng)式對(duì)這一組數(shù)據(jù)適用，對(duì)另一組數(shù)據(jù)可能就要用另外的多項(xiàng)式.因此，在實(shí)際問題中使用哪種方法要視具體問題而定，有時(shí)幾種方法都要用到.實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分小結(jié)積分描述了一個(gè)函數(shù)的整體或者宏觀性質(zhì)，對(duì)函數(shù)的形狀在小范圍的改變不敏感，并且積分過程是對(duì)數(shù)據(jù)點(diǎn)進(jìn)行求和，正的和負(fù)的隨機(jī)誤差傾向于相互抵消，因此，一般數(shù)值積分過程是穩(wěn)定的，所得的解精確度也較高.微分則描述了一個(gè)函數(shù)在一點(diǎn)處的斜率，是函數(shù)的微觀性質(zhì)，它很敏感，一個(gè)函數(shù)小的變化，容易產(chǎn)生相鄰點(diǎn)的斜率的大的改變.8586實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分差商近似就是用給定函數(shù)f(x)曲線上的割線斜率近似切線斜率，插值型數(shù)值微分和擬合型數(shù)值微分都是用近似多項(xiàng)式的曲線斜率近似給定函數(shù)f(x)的曲線斜率.無論是割線斜率，還是近似多項(xiàng)式的曲線斜率，都可能和給定函數(shù)f(x)的曲線斜率有很大不同，特別是當(dāng)f(x)在給定區(qū)間內(nèi)變化比較大時(shí)更是這樣.同時(shí)，微分過程是對(duì)數(shù)據(jù)點(diǎn)進(jìn)行相減，正的和負(fù)的隨機(jī)誤差傾向于相加，這都使得數(shù)值微分的解不穩(wěn)定，并且精度也較差.第四章數(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ù)值解88實(shí)驗(yàn)4.4常微分方程的數(shù)值解一、幾種求常微分方程數(shù)值解的方法二、應(yīng)用舉例實(shí)驗(yàn)4.4常微分方程的數(shù)值解89實(shí)驗(yàn)4.4常微分方程的數(shù)值解實(shí)驗(yàn)?zāi)康耐ㄟ^本實(shí)驗(yàn)了解常微分方程的數(shù)值解的概念，掌握利用MATLAB求常微分方程的數(shù)值解的方法.一、幾種求常微分方程數(shù)值解的方法常微分方程是研究函數(shù)變化規(guī)律的有力工具，但是絕大多數(shù)變系數(shù)方程、非線性方程都是所謂“解不出來”的，于是常微分方程的數(shù)值解法就成為解常微分方程的主要手段．本實(shí)驗(yàn)只考慮初值問題.90實(shí)驗(yàn)4.4常微分方程的數(shù)值解常微分方程初值問題：保證方程的解存在且惟一在一系列離散點(diǎn)上，求的近似值，通常取等步長(zhǎng)h，即因此數(shù)值解法得到的近似解是一個(gè)離散的函數(shù)表.91實(shí)驗(yàn)4.4常微分方程的數(shù)值解一、幾種求常微分方程數(shù)值解的方法1.歐拉（Euler）方法歐拉方法是一種最古老、最簡(jiǎn)單而直觀的解微分方程的數(shù)值方法，其基本想法是在小區(qū)間上用差商代替方程左端的導(dǎo)數(shù),而方程右端已知函數(shù)中的x在小區(qū)間上的哪一點(diǎn)取值，則有以下不同的方法：92實(shí)驗(yàn)4.4常微分方程的數(shù)值解（1）向前歐拉公式中的x取小區(qū)間的左端點(diǎn)xn，以分別代替，就得到：向前歐拉公式（2）向后歐拉公式中的x取小區(qū)間的右端點(diǎn)xn+1，就得到：向后歐拉公式上式右端的yn+1未知，故稱為隱式公式，無法用它直接計(jì)算yn+1．93實(shí)驗(yàn)4.4常微分方程的數(shù)值解（3）梯形公式將向前歐拉公式和向后歐拉公式加以平均，得到梯形公式94實(shí)驗(yàn)4.4常微分方程的數(shù)值解（4）改進(jìn)的歐拉公式先由向前歐拉公式算出yn+1的預(yù)測(cè)值，再把它代入梯形公式右端，作為校正，即改進(jìn)的歐拉公式95實(shí)驗(yàn)4.4常微分方程的數(shù)值解它還可寫作上面4個(gè)公式中，我們常用的是便于計(jì)算的向前歐拉公式和改進(jìn)的歐拉公式.96實(shí)驗(yàn)4.4常微分方程的數(shù)值解2.龍格—庫(kù)塔（Runge-Kutta）公式在區(qū)間內(nèi)多取幾個(gè)點(diǎn)，就可以構(gòu)造出精度更高的計(jì)算公式，這就是龍格—庫(kù)塔公式，它是一類方法的總稱.它是歐拉方法的一種推廣，也是應(yīng)用最廣的求解常微分方程數(shù)值問題的方法.一般的龍格—庫(kù)塔方法的形式為p階龍格—庫(kù)塔方法其中ai,bij,ci為待定參數(shù)97實(shí)驗(yàn)4.4常微分方程的數(shù)值解（1）二階龍格—庫(kù)塔公式其中為待定系數(shù).滿足：上式有4個(gè)未知數(shù)而只有3個(gè)方程，所以解不唯一，存在無窮多個(gè)解，可見二階龍格—庫(kù)塔方法是一族公式.

改進(jìn)的歐拉公式98實(shí)驗(yàn)4.4常微分方程的數(shù)值解（2）四階經(jīng)典（標(biāo)準(zhǔn)）龍格—庫(kù)塔公式實(shí)際應(yīng)用中，用的最多的是四階龍格—庫(kù)塔公式四階龍格—庫(kù)塔公式也不止一個(gè)99實(shí)驗(yàn)4.4常微分方程的數(shù)值解3.龍格—庫(kù)塔方法的MATLAB實(shí)現(xiàn)在MATLAB中，求微分方程數(shù)值解的函數(shù)調(diào)用格式為[t,x]=solver（’f’,ts,x0,options）自變量值函數(shù)值ode45、ode23、ode113、ode15s、ode23s中的一個(gè)由待解方程寫成的m-文件名=[t0,tf],t0、tf為自變量的初值和終值函數(shù)的初值用于設(shè)定誤差限options=odeset（'reltol',rt,'abstol',at）默認(rèn)相對(duì)誤差為10-3，絕對(duì)誤差為10-6.相對(duì)誤差絕對(duì)誤差100實(shí)驗(yàn)4.4常微分方程的數(shù)值解在諸多的solver函數(shù)中，ode23、ode45都是基于龍格—庫(kù)塔公式的函數(shù):ode23是采用的2、3階龍格-庫(kù)塔組合算法，ode45是采用的4、5階龍格-庫(kù)塔組合算法.ode45比ode23精度更高，尤其成為大部分場(chǎng)合的首選算法.注(1)在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成.(2)使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.101例21

求微分方程的數(shù)值解，解先編寫函數(shù)文件：functiondx=fun23(t,x)dx=t+x;再求方程的解實(shí)驗(yàn)4.4常微分方程的數(shù)值解已知精確解為>>ts=0:0.1:1;x0=1;[t,x]=ode45('fun23',ts,x0);y=2*exp(t)-t-1;[t,x,y]plot(t,x,'r-',t,y,'b-.'),gridon;title('SolutionofExample3');xlabel('timet');ylabel('solutionx');legend('x','y');↙ans=01.00001.0000102實(shí)驗(yàn)4.4常微分方程的數(shù)值解0.10001.11031.11030.20001.24281.24280.30001.39971.39970.40001.58361.58360.50001.79741.79740.60002.04422.04420.70002.32752.32750.80002.65112.65110.90003.01923.01921.00003.43663.4366

t的值x的值y的值例22

求微分方程組的數(shù)值解.實(shí)驗(yàn)4.4常微分方程的數(shù)值解103解先編寫函數(shù)文件：functiondx=fu