第六章數(shù)值計算命令與例題

1、 第六章 數(shù)值計算命令與例題數(shù)值計算命令與例題l北京交通大學(xué)6.1 求近似函數(shù)求近似函數(shù)l在生產(chǎn)和實驗中在生產(chǎn)和實驗中, 人們經(jīng)常遇到需要通過某個未知的函數(shù)人們經(jīng)常遇到需要通過某個未知的函數(shù)f(x)在有限個給定點的函數(shù)值在有限個給定點的函數(shù)值:xi, yi, i=1,2,., n, 這里這里 f(xi) = yi 去獲得函數(shù)去獲得函數(shù)f(x)的近似函數(shù)的近似函數(shù) (x), 求近似函數(shù)求近似函數(shù) (x)的方法的方法主要有主要有擬合擬合方法和方法和插值插值方法。方法。6.1.1 曲線擬合曲線擬合l 曲線擬合主要用來求一元近似函數(shù), 它是根據(jù)最小二乘原理的意義下獲得近似函數(shù)的, 此近似函數(shù)具有在數(shù)據(jù)

2、點處的誤差平方和最小的特點。記函數(shù)集合:M=Span 0, 1, 2, m= (x)| (x)= a0 0(x)+a1 1(x)+am m(x), ai Rl稱集合稱集合M為函數(shù)為函數(shù) 0, 1, 2, m張成的空間，張成的空間，m+1個函數(shù)個函數(shù) 0(x), 1(x), 2(x), m(x)稱為擬合基函數(shù)集合稱為擬合基函數(shù)集合, 它們都是它們都是已知的函數(shù)。已知的函數(shù)。lMathematica曲線擬合的一般形式為曲線擬合的一般形式為: Fit數(shù)據(jù)點集合數(shù)據(jù)點集合, 擬合基函數(shù)集合擬合基函數(shù)集合, 自變量名自變量名具體的擬合命令有具體的擬合命令有: 命令形式命令形式1：Fitx1,y1,x2,

3、y2,.,xn,yn, 0, 1, 2, m ,x功能：根據(jù)數(shù)據(jù)點集功能：根據(jù)數(shù)據(jù)點集x1,y1,x2,y2,.,xn,yn求出具有擬合函數(shù)為求出具有擬合函數(shù)為 (x)= a 0 0(x)+a1 1(x)+a m m(x) 形式的近似函數(shù)形式的近似函數(shù) (x) 命令形式命令形式2：Fity1,y2,.,yn, 0, 1, 2, m ,x功能：根據(jù)數(shù)據(jù)點集功能：根據(jù)數(shù)據(jù)點集1,y1,2,y2,.,n,yn求出具有擬合函數(shù)為求出具有擬合函數(shù)為 (x)= a 0 0(x)+a1 1(x)+a m m(x)形式的近似函數(shù)形式的近似函數(shù) (x)命令形式命令形式3：Fitx1,y1,x2,y2,.,xn,

4、 yn, Tablexi,i,0,m ,x功能：根據(jù)數(shù)據(jù)點集功能：根據(jù)數(shù)據(jù)點集x1,y1,x2,y2,.,xn,yn求出擬合函數(shù)為求出擬合函數(shù)為m次多次多項式的近似函數(shù)項式的近似函數(shù) (x) =a 0 +a1x+ a2x2 +a mx m多項式擬合算法多項式擬合算法 l輸入n+1個擬合點： (xi, yi),i=0,1,nl根據(jù)散點圖確定擬合多項式的次數(shù)ml計算相應(yīng)正規(guī)線性方程組的系數(shù)和右端項l解正規(guī)正規(guī)線性方程組，得解：a0*,a1*,a m*l寫出擬合多項式*(x)= a0*+ a1*x+ a2*x2+ + am*xm求求m次多項式擬合程序次多項式擬合程序lClearxi,xx,yi;lx

5、i=Inputxi=lyi=Inputyi=ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04lm=Input多項式次數(shù)m=ls=TableSumxiki,k,1,n,i,0,2m;la=Tablesi+j-1,i,1,m+1,j,1,m+1;lPrinta=,MatrixForma;lb=TableSumxiki*yik,k,1,n,i,0,m;lPrintb=,MatrixFormb;lxx=Tablexi,i,1,m+1;lg=Solvea.xx=b,xx;lfa=Sumxi*t(i-1),i,1,m+1/.

6、g1;lp=fa/Nlp1=Plotp,t,xi1,xin,DisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunction;l說明：說明：本程序用于求m次多項式擬合。程序執(zhí)行后，按要求通過鍵盤輸入擬合基點xi:x0 , x1, . , xn 、對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，計算機給出散點圖和請求輸入擬合多項式次數(shù)的窗口，操作者可以根據(jù)散點圖確定擬合多項式的次數(shù)通過鍵盤輸入，程序即可給出對應(yīng)的正規(guī)方程組系數(shù)矩陣a、常數(shù)項b、m次擬合多項式和由擬合函數(shù)圖形和散點圖畫在一起的圖形。l程序中變量說明程序中變量

7、說明lxi：存放擬合基點x0 , x1, . , xn lyi: 存放對應(yīng)函數(shù)值y0 , y1 , , ynlm: 存放擬合多項式次數(shù)la: 存放正規(guī)方程組系數(shù)矩陣lb: 存放正規(guī)方程組常數(shù)項lp: 存放m次擬合多項式lh: 存放散點圖lp1：存放擬合函數(shù)圖形lxx:定義正規(guī)方程組變量，存放m次擬合多項式的系數(shù)l注：注：語句s=TableSumxiki,k,1,n,i,0,2m、a=Tablesi+j-1,i,1,m+1,j,1,m+1、lb=TableSumxiki*yik,k,1,n,i,0,m是用簡化的正規(guī)方程組編程的。例例1已知一組實驗數(shù)據(jù)x 1 3 4 5 6 7 8 9 10f(x

8、) 10 5 4 2 1 1 2 3 4用多項式擬合求其擬合曲線。l解：解：執(zhí)行m次多項式擬合程序后，在輸入的兩個窗口中按提示分別輸入1,3,4,5,6,7,8,9,10，10,5,4,2,1,1,2,3,4l每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，計算機在屏幕上畫出散點圖。l由于該散點圖具有2次多項式形狀，因此在確定選擇多項式次數(shù)窗口輸入2，按OK”按扭后得如下輸出結(jié)果。la= 9 53 381 53 381 3017 381 3017 25317lb=32147102513.4597 - 3.60531 t + 0.267571 t2l于是得求出的二次多項式擬合函數(shù)為(t)=13.4597

9、 - 3.60531 t + 0.267571 t2而且從圖形上看擬合效果很好。l例例 2已知一組實驗數(shù)據(jù)lx 6 8 10 12 14 16 18 20 22 24lf(x) 4.6 4.8 4.6 4.9 5.0 5.4 5.1 5.5 5.6 6.0l修改多項式擬合程序使其具有可以選擇不同多項式擬合函數(shù)功能，并用此程序確定本題多項式擬合曲線。l解：解： 在擬合多項式程序中加入評價擬合效果的判定人機交互語句 tu=InputSatisfatory?0(No)or1Yes 和While語句來調(diào)控何時終止實驗，調(diào)控變量用tu取值確定，取值0表示不滿意，1表示滿意。此外，去掉正規(guī)方程組系數(shù)及擬合

10、多項式的輸出，代之以在圖形上表出擬合多項式的次數(shù)提示，修改后的程序如下。lxi=Tablei,i,6,24,2;lyi=4.6,4.8,4.6,4.9,5,5.4,5.1,5.5,5.6,6;ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04ltu=0;lWhiletu=0,lm=Input多項式次數(shù)m;ls=TableSumxiki,k,1,n,i,0,2m;la=Tablesi+j-1,i,1,m+1,j,1,m+1;l(*Printa=,MatrixForma;*)lb=TableSumxiki*yik,k,

11、1,n,i,0,m;l(*Printb=,MatrixFormb;*)lxx=Tablexi,i,1,m+1;lg=Solvea.xx=b,xx;lfa=Sumxi*t(i-1),i,1,m+1/.g1;lp=fa/N;lp1=Plotp,t,xi1,xin,PlotLabel-m擬合多項式,lDisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunction;ltu=InputSatisfatory?0(No)or1Yesl執(zhí)行該程序后，屏幕上出現(xiàn)擬合多項式次數(shù)輸入窗口和散點圖。l由于該散點圖不好確定擬合次數(shù)，先用3次擬合多

12、項式計算，因此輸入“3”并用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出圖形。l屏幕出現(xiàn)提示是否滿意的輸入窗口，因為不太滿意，輸入：0，單擊“OK”按扭并在屏幕上出現(xiàn)擬合多項式次數(shù)輸入窗口中輸入：6，單擊OK，屏幕出現(xiàn)下一個組合圖形。l繼續(xù)做實驗，得到如下若干圖形。l由于總共有10個數(shù)據(jù)點，所以擬合多項式最高次數(shù)只能到9次，因此實驗到9次擬合多項式后，在屏幕出現(xiàn)提示是否滿意的輸入窗口中，輸入：1，單擊“OK”按扭退出實驗。線性模型擬合線性模型擬合 l線性模型擬合算法線性模型擬合算法l l1.輸入n+1個擬合點： (xi, yi),i=0,1,nl2.根據(jù)散點圖確定擬合基函數(shù)組l3.計算相應(yīng)正規(guī)線性方

13、程組的系數(shù)和右端項l4.解正規(guī)正規(guī)線性方程組，得解：a0*,a1*,a m*l5.寫出線性擬合模型*(x)= a0*0(x)+ a1*1(x)+ a2*2(x)+ + am*m(x)求線性模型擬合程序求線性模型擬合程序lClearx,xi,xx,yi;lxi=Inputxi=lyi=Inputyi=ln=Lengthxi;lh=ListPlotTablexii,yii,i,1,n,PlotStyle-PointSize0.04lm1=Input輸入擬合基函數(shù)組lm=Lengthm1lp=Tablem1/.x-xik,k,1,nla=TableSumpk,i*pk,j,k,1,n,i,1,m,j

14、,1,m/Nl(*Printa=,MatrixForma;*)lb=TableSumpk,i*yik,k,1,n,i,1,m/Nl(*Printb=,MatrixFormb;*)lxx=Tablexti,i,1,mlg=Solvea.xx=b,xxlfa=Sumxti*m1i,i,1,m/.g1lpp=fa/N;lp1=Plotpp,x,xi1,xin,DisplayFunction-Identity;lShowp1,h,DisplayFunction-$DisplayFunctionl說明：說明：本程序用于求線性模型擬合。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi:x0 , x1, . ,

15、xn 、對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，計算機給出散點圖和請求輸入擬合擬合基函數(shù)組0(x)，1(x)，2(x)、m(x)的窗口，操作者可以根據(jù)散點圖確定擬合基函數(shù)組通過鍵盤輸入，程序即可給出對應(yīng)的正規(guī)方程組系數(shù)矩陣a、常數(shù)項b、線性模型擬合函數(shù)和由擬合函數(shù)圖形與散點圖畫在一起的圖形。l程序中變量說明：程序中變量說明：lxi：存放擬合基點x0 , x1, . , xn lyi: 存放對應(yīng)函數(shù)值y0 , y1 , , ynlm1: 存放擬合基函數(shù)組0(x)，1(x)，2(x)、m(x)lm: 存放擬合基函數(shù)組個數(shù)la: 存放正規(guī)方程組系數(shù)矩陣lb: 存放正規(guī)方程組常數(shù)項lp:

16、存放擬合基函數(shù)組在擬合基點x0 , x1, . , xn 的函數(shù)值lpp: 存放求出的線性模型擬合函數(shù)lh: 存放散點圖lp1：存放擬合函數(shù)圖形lxx:定義正規(guī)方程組變量，存放線性模型擬合系數(shù)l注：（注：（1）語句m1=Input輸入擬合基函數(shù)組產(chǎn)生的輸入應(yīng)用函數(shù)表，即用l“0 x，1x，mx”的形式輸入。l（2）Mathematica數(shù)學(xué)軟件中也有一個求線性模型擬合的命令，命令形式為l Fitx1,y1,x2,y2,.,xn,yn, 0, 1, 2, m ,xl它可以根據(jù)數(shù)據(jù)點集x1,y1,x2,y2,.,xn,yn求出具有擬合函數(shù)為l (x)= a 0 0(x)+a11(x)+a mm(x

17、)l形式的近似函數(shù)(x)。l例例3已知函數(shù)在自變量x=1,2, , 10上數(shù)據(jù)為2.89229, 2.86323, 0.473147, -2.25209, -2.87003, -0.835768, 1.97187, 2.96841, 1.23648, -1.63202，試用合適的函數(shù)進行擬合。l解：解：執(zhí)行線性模型擬合程序后，在輸入的兩個窗口中按提示分別輸入l1,2,3,4,5,6,7,8,9,10、2.89229, 2.86323, 0.473147, -2.25209, -2.87003, -0.835768, 1.97187, 2.96841, 1.23648, -1.63202l每次輸

18、入后用鼠標(biāo)點擊窗口的“OK”按扭，計算機在屏幕上畫出散點圖。l由于該散點圖具有類似正弦曲線的形狀，因此在確定選擇擬合基函數(shù)窗口輸入“1,Sinx”，按“OK”按扭后得如下輸出結(jié)果。la=10. 1.41119l 1.41119 5.00143lb=4.81552l 15.4239l0.0482789 + 3.07027 Sinxl于是，我們得到了很好的擬合函數(shù)(x)=0.0482789 + 3.07027 sin(x) 。l對于本題，如果看到散點圖具有一個彎曲而選用三次多項式擬合，則輸入“1,x,x2,x3”后會得到如下輸出。la=l 10. 55. 385. 3025.l 55. 385.

19、3025. 25333.l 385. 3025. 25333. 220825.l 3025. 25333. 220825. 1.97841 106lb=4.81552l 14.0236l 104.284l 820.711l10.4765 - 7.37686 x + 1.41989 x2 - 0.0796299 x3l顯然這個擬合很不滿意。6.1.2 函數(shù)插值函數(shù)插值l多項式插值多項式插值 多項式插值是在給定多項式插值是在給定n個數(shù)據(jù)點個數(shù)據(jù)點:xi, yi, i=1,2,., n 后后, 求出一個次求出一個次數(shù)不超過數(shù)不超過n-1的多項式的多項式 (x)作為函數(shù)作為函數(shù) y=f(x) 的近似表

20、達(dá)式，它就是的近似表達(dá)式，它就是計算方法中常說的計算方法中常說的拉格朗日拉格朗日(Lagrange)插值或插值或Newton插值多項式。插值多項式。lMathematica多項式插值的一般形式為多項式插值的一般形式為: InterpolatingPolynomial 數(shù)據(jù)點集合數(shù)據(jù)點集合, 自變量名自變量名l具體的多項式插值命令有具體的多項式插值命令有:命令形式命令形式1：InterpolatingPolynomial x1,y1,x2,y2,.,xn,yn,x功能：功能：根據(jù)數(shù)據(jù)點集x1,y1,x2,y2,.,xn,yn求出n-1次插值多項式 (x)命令形式命令形式2：Interpolati

21、ngPolynomial y1,y2,.,yn,x功能：功能：根據(jù)數(shù)據(jù)點集1,y1,2,y2,.,n,yn求出n-1次插值多項式 (x)Lagrange插值插值l求求n次次Lagrange插值多項式算法插值多項式算法l l輸入n+1個插值點： (xi, yi),i=0,1,nl計算插值基函數(shù)l0n(x), l1n(x) , lnn(x)l給出n次Lagrange插值多項式：Ln(x)= y0 l0n(x)+ y1 l1n(x) +yn lnn(x)l求求Lagrange插值多項式程序插值多項式程序lClearlag,xi,x,yi;lxi=Inputxi=lyi=Inputyi=ln=Leng

22、thxi-1;lp=Sumyii*(Product(x-xij)/(xii-xij),j,1,i-1l*Product(x-xij)/(xii-xij),j,i+1,n+1),i,1,n+1;llagx_=Simplifypl說明：說明：本程序用于求n次Lagrange插值多項式。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi:x0 , x1, . , xn 和對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，程序即可給出對應(yīng)的n次Lagrange插值多項式lagx。l程序中變量說明程序中變量說明lxi：存放插值基點x0 , x1, . , xn lyi: 存放對應(yīng)函數(shù)值y0 , y1 , ,

23、ynllagx: 存放求出的n次Lagrange插值多項式Ln(x)l注：注：語句lagx_=Simplifyp用簡化形式給出對應(yīng)的n次Lagrange插值多項式。l例例給定數(shù)據(jù)表 x 0 1 2 3 y=f(x) 1 3 5 12 用Lagrange插值法求三次插值多項式,并給出函數(shù)f(x)在x =1.4的近似值。l解：解： 執(zhí)行Lagrange插值程序后，在輸入的兩個窗口中按提示分別輸入0, 1, 2, 3、1, 3, 5, 12，每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)。 6 + 22 x - 15 x2 + 5 x3 - 6l所以得到三次插值多項式L3(x)=1+11 x

24、/3-5 x2/2+5 x3/6 接著鍵入“l(fā)ag1.4”，則輸出3.52，因此f(x)在x =1.4的近似值為3.52，即f(1.4)3.52.Newton插值插值l 求求Newton插值多項式算法插值多項式算法l l1. 輸入n+1個插值點： (xi, yi),i=0,1,nl2. 計算差商表l3.給出n次Newton插值多項式。求求Newton插值多項式程序插值多項式程序lClearnewt,s,x;lxi=Inputxi=lyi=Inputyi=ln=Lengthxi;l(*計算差商表*)lf=Table0,n,n;lDofi,1=yii,i,1,nlDofi,j+1=(fi,j-fi

25、+1,j)/(xii-xii+j), j,1,n-1,i,1,n-jlPrint差商表lDoPrintxii, ,fi,i,1,nl(*求Newton插值多項式*)lfa=1;ls=f1,1;lDofa=(x-xik)*fa;s=s+fa*f1,k+1,k,1,n-1lnewtx_=slSimplify%l說明：說明：本程序用于求n次Newton插值多項式。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi:x0 , x1, . , xn 和對應(yīng)函數(shù)值yi: y0 , y1 , , yn 后，程序依次給出輸入的數(shù)據(jù)表、計算出的差商表、Newton插值多項式、Newton插值多項式的簡化形式。l程序中變

26、量說明：程序中變量說明：lxi：存放插值基點x0 , x1, . , xn lyi: 存放對應(yīng)函數(shù)值y0 , y1 , , ynlf:存放函數(shù)值y0 , y1 , , yn及所有差商lnewtx: 存放求出的n次newton插值多項式Nn(x)l注：注：l（1)語句f=Table0,n,n用于產(chǎn)生一個nn的矩陣變元用于存放函數(shù)值y0 , y1 , , yn及所有差商。l（2)在Mathematica中有一個求n次插值多項式的命令，命令形式 InterpolatingPolynomialx01,y0，x1,y1,x2,y2,xn,yn, xl它可以求過n+1個插值點x01,y0，x1,y1,x2

27、,y2,xn,yn的n次插值多項式Pn(x)。l例例2給定數(shù)據(jù)表l x 4.0002 4.0104 4.0233 4.0294l y 0.6020817 0.6031877 0.6045824 0.6052404l(1)計算差商表l(2)用Newton插值法求三次插值多項式Nn(x)l(3)求f(4.011)的近似值l解：解： 執(zhí)行Newton插值程序后，在輸入的兩個窗口中按提示分別輸入4.0002, 4.0104, 4.0233, 4.0294，0.6020817, 0.6031877, 0.6045824, 0.6052404,每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)。4.0

28、002, 4.0104, 4.0233, 4.02940.6020817, 0.6031877, 0.6045824, 0.6052404 l差商表4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.0130225, 04.0233 0.6045824, 0.107869, 0, 04.0294 0.6052404, 0, 0, 00.6020817 + 0.108431 (-4.0002 + x) - 0.0136404 (-4.0104 + x) (-4.0002 + x) + 0.

29、0211629 (-4.0233 + x) (-4.0104 + x) (-4.0002 + x)-1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3l于是我們得到本題的差商表為xi yi f, f , , f , , ,4.0002 0.6020817, 0.108431, -0.0136404, 0.02116294.0104 0.6031877, 0.108116, -0.01302254.0233 0.6045824, 0.107869 4.0294 0.6052404 和三次插值多項式N3(x)= -1.41642 + 1.23926 x

30、 - 0.268313 x2 + 0.0211629 x3接著鍵入“newt4.011”，則輸出0.603253，因此f(4.011) 0.603253。l例例3多項式插值的誤差估計式中可以看到，當(dāng)插值節(jié)點越多時誤差會越小，這個結(jié)論正確嗎？通過實驗說明該結(jié)論的正確性。l解解: 考慮函數(shù)f(x) = (1+x2)-1 在區(qū)間-4,4內(nèi)選取不同個數(shù)的等距插值節(jié)點做觀察，這里分別選-4,4 內(nèi)的9個和11個的等距節(jié)點來做實驗，將對應(yīng)的插值函數(shù)圖與被插函數(shù)f(x) = (1+x2)-1畫在一起做觀察，為簡單起見，這里用Mathematica 命令做實驗，對應(yīng)命令為lIn1:= u=Tablex,(1+

31、x2)-1,x,-4,4 ; （*采取f(x) 在-4,4 內(nèi)的9個插值點 lIn2:= g=ListPlotu, PlotStyle-PointSize0.04 （*將散點圖圖形文件存放在變量g中l(wèi)In3:= s=InterpolatingPolynomialu , x ; （*將插值函數(shù)存放在變量s中l(wèi)In4:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle-Thickness0.005,Thickness0.006 （*將插值函數(shù)s與f(x)畫在一起的圖形文件存放在變量t中l(wèi)In5:= Showt, g（*將散點圖, 插值函數(shù)s, f(x)畫在一個坐標(biāo)系中l(wèi)

32、在-4,4中選9個等距節(jié)點的插值函數(shù)與被插函數(shù)圖，粗線為被插函數(shù)圖 lIn6:= u=Tablex,(1+x2)-1,x,-4,4，0.8 ; （*采取f(x) 在-4,4 內(nèi)的11插值點 lIn7:= g=ListPlotu, PlotStyle-PointSize0.04 lIn8:= s=InterpolatingPolynomialu , x ; lIn9:= t= Plots, (1+x2)-1, x,-4,4, PlotStyle-Thickness0.005,Thickness0.006 lIn10:= Showt, g l在-4,4中選11個等距節(jié)點的插值函數(shù)與被插函數(shù)圖 分段

33、多項式插值分段多項式插值Mathematica分段多項式插值的一般形式為分段多項式插值的一般形式為: Interpolation 數(shù)據(jù)點集合數(shù)據(jù)點集合具體的分段多項式插值命令有具體的分段多項式插值命令有:l命令形式命令形式1：Interpolation x1,y1,x2,y2,.,xn,yn功能：功能：根據(jù)數(shù)據(jù)點集x1,y1,x2,y2,.,xn,yn求出分段插值多項式(x)l命令形式命令形式2：Interpolation y1,y2,.,yn功能：功能：根據(jù)數(shù)據(jù)點集1,y1,2,y2,.,n,yn求出分段插值多項式(x)例5: 用分段多項式插值重做例4的插值問題, 并驗證, 隨著插值點的增多

34、, 分段插值函數(shù)的誤差會越來越小。(見圖) 解解:In27:= r=InterpolationTablex,(1+x2)-1, x,-4,4 Out27= InterpolatingFunction-4, 4, In28:= (1+x2)-1/.x-3.7Out28= 0.0680735 In29:= r3.7Out29=0.0734 In30:= d= Tablex,(1+x2)-1, x,-4,4,0.5; In31 := r1=InterpolationdOut31= InterpolatingFunction-4, 4., In32:= r13.7Out32= 0.0681761In3

35、3:= Plotr1x, (1+x2)-1, x,-4,4 下一頁下一頁返回返回分段線性插值分段線性插值l用分段線性插值函數(shù)做近似計算的算法用分段線性插值函數(shù)做近似計算的算法l1. 輸入n個插值點： (xi, yi),i=1,nl2. 輸入要近似計算的自變量點xal3.尋找包含xa的小區(qū)間xi , xi+1l4.用 xi , xi+1 上的線性插值函數(shù)在xa處的值作為f(xa)的函數(shù)值l用分段線性插值函數(shù)做近似計算的程序用分段線性插值函數(shù)做近似計算的程序lClearx,a,b;llia_,b_,x_:=(b2-a2)/(b1-a1)*(x-a1)+a2lxi=Inputxi=lyi=Input

36、yi=lxa=Inputxa=ln=Lengthxi;lIfxaxin,Print超限;Break,l DoIfxa=xik,m=k,k,1,n-1;l Printm=,m,x=,xim;l lixim,yim,xim+1,yim+1,xal說明：說明：本程序用分段線性插值做近似計算。程序執(zhí)行后，按要求通過鍵盤輸入插值基點xi: x1, . , xn 和對應(yīng)函數(shù)值yi: y1 , , yn 和要做近似計算的點xa后，程序?qū)⒔o出包含xa小區(qū)間xi , xi+1的下標(biāo)i和xi ,和函數(shù)f(xa)的近似值。如果xa不在x1, xn 內(nèi)程序給出超限提示。l程序中變量說明程序中變量說明lxi：存放插值基

37、點 x1, . , xn lyi: 存放對應(yīng)函數(shù)值 y1 , , ynlxa: 存放要做近似計算的自變量值lm:存放包含xa小區(qū)間xi , xi+1的下標(biāo)il注：在Mathematica中有一個求分段線性插值函數(shù)的命令，命令形式為l Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder - 1l用如上Mathematica命令求出的分段線性插值函數(shù)沒有給出具體的分段函數(shù)表達(dá)式, 而是用l“InterpolatingFunctionx1, xn, ”作為所求的分段插值函數(shù), 通?？梢杂胠 變量=Interpolation 數(shù)據(jù)點集合l把所得的分

38、段插值函數(shù)存放在變量中, 如果要計算插值函數(shù)在某一點的如p點的值,只要輸入 “變量p” 即可，如果要表示這個插值函數(shù)應(yīng)該用 “變量x” 。l此外，命令l Interpolation x1,y1,x2,y2,.,xn,yn，InterpolationOrder - 3l可以給出更好的分段插值函數(shù)。l例例4給定數(shù)據(jù)表 xi 1 2 3 4 yi 0 -5 -6 3l（1）用分段線性插值函數(shù)求函數(shù)f(x)在x =1.4的近似值l（2）用Mathematica命令畫出分段線性插值圖形l解：解：執(zhí)行求分段線性插值函數(shù)程序后，在輸入的兩個窗口中按提示分別輸入l1，2，3，4，0，-5，-6，3，1.4l每

39、次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下插值函數(shù)。lm=1 x=1l-2l此結(jié)果說明f(x)在x =1.4的近似值為2，它對應(yīng)的小區(qū)間為x1 , x2。為完成第二個問題，鍵入命令ld=Interpolation 1,0,2,-5,3,-6,4,3,InterpolationOrder - 1lPlotdx,x,1,4l得輸出如下。lInterpolatingFunction1, 4, 6.2 解常微分方程l自變量、未知函數(shù)以及未知函數(shù)的導(dǎo)數(shù)（或微分）之間的一個（或一組）關(guān)系式稱為微分方程（或微分方程組）。求出常微分方程（或微分方程組）的未知函數(shù)（或未知函數(shù)組），稱為解常微分方程。含有任意常

40、數(shù)的解稱為通解，否則稱為特解。n階常微分方程的一般形式為：l利用Mathematica可以象高等數(shù)學(xué)一樣求出常微分方程的解析解（準(zhǔn)確解），如果求不出解析解Mathematica可以求出微分方程的數(shù)值解（即近似解。Mathematica解常微分方程的命令有求常微分方程（組）的解析解命令和求常微分方程（組）的數(shù)值解命令。 0),()(nyyyxF6.2.1求常微分方程（組）的解析解求常微分方程（組）的解析解命令形式命令形式1：DSolveeqn, yx, x 功能：功能：求出常微分方程eqn的未知函數(shù)yx的解析通解。命令形式命令形式2：DSolveeqn1,eqn2, . , y1x,y2x, .

41、 , x 功能：功能：求出常微分方程組eqn1,eqn2, . 的所有未知函數(shù)y1x,y2x, . 的解析通解。注意注意：l每個常微分方程的中的等號輸入時應(yīng)該用兩個等號代替l未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)要用“自變量名”指示未知函數(shù)的自變量l常微分方程組和未知函數(shù)組應(yīng)該用大括號括起來l命令形式2可以用來求常微分方程的特解leqn表示常微分方程, x是自變量名, yx 表示未知函數(shù), 自變量名和未知函數(shù)可以是其他的變量名例例6: 求常微分方程y =2ax的通解，a為常數(shù)解:Mathematica 命令: In34:= DSolveyx = 2 a x, yx, xOut34= yx - a x 2 +

42、 C1即得本問題的通解 例7: 求常微分方程y +y =0的通解 解:Mathematica 命令: In35:= DSolveyx+yx=0, yx, xOut35= yx - C2 Cosx - C1 Sinx即得本問題的通解:C1，C2是任意常數(shù)。 xcxcysincos1221xacy例8: 求常微分方程的特解。 解:Mathematica 命令: In36:= DSolveyx=x/yx+yx/x, y1=2, yx, x Out36= yx - Sqrtx2 (4 + 2 Logx) 本問題的特解為:l下面的操作可以檢驗所求特解的正確性:In37:= yx_:= Sqrtx2* (

43、4 + 2 Logx) In38:= y1 Out38=2In39:=Dyx,x-x/yx-yx/x; In40:= Simplify% Out40= 0 2|,1xyxyyxy)ln24(2xxy例例9: 求常微分方程組求常微分方程組y (t) =x(t), x (t) =y(t) 的通解。的通解。 解解:Mathematica 命令命令: In41:= DSolveyt=xt,xt=yt, xt,yt,t C1+ E 2 t C1-C2+ E 2 t C2l Out41= xt - -, 2 E t -C1+ E 2 t C1+C2+ E 2 t C2 yt - - 2 E t6.2.2

44、求常微分方程（組）的數(shù)值解命令求常微分方程（組）的數(shù)值解命令l命令形式命令形式1：NDSolve eqns, y, x, xmin, xmax 功能：功能：求出自變量范圍為xmin, xmax且滿足給定常微分方程及初值條件eqns的未知函數(shù)y的數(shù)值解 l命令形式命令形式2：NDSolve eqns, y1,y2, . , x, xmin, xmax 功能：功能：求出自變量范圍為xmin, xmax且滿足給定常微分方程及初值條件eqns的未知函數(shù)y1,y2, . 的數(shù)值解。例例10:求微分方程初值問題求微分方程初值問題y =x+y,y(0)=1在區(qū)間在區(qū)間0,0.5內(nèi)的數(shù)值解內(nèi)的數(shù)值解解解:Ma

45、thematica 命令命令: In42:= NDSolve yx=x+yx, y0=1, y, x,0,0.5 Out42= y - InterpolatingFunction0., 0.5, l上面顯示的是所求數(shù)值解的替換形式，為得到本問題數(shù)值解，再鍵入上面顯示的是所求數(shù)值解的替換形式，為得到本問題數(shù)值解，再鍵入:In43:= f= y/.% 1 Out43= InterpolatingFunction0., 0.5, In44:= Plotfx,x,0,0.5 In45:= Table fx, x, 0, 0.5, 0.1 Out45= 1., 1.11034, 1.24281, 1.3

46、9972, 1.58365, 1.79744 例例11: 已知常微分方程組已知常微分方程組求函數(shù)求函數(shù)x(t)和和y(t)在在0,10上的數(shù)值解。上的數(shù)值解。解解:Mathematica 命令命令: In46:=q=NDSolvext=-xt2-yt,yt=2xt-yt,x0=y0=1,x,y,t,0,10 Out46= x - InterpolatingFunction0., 10., , y - InterpolatingFunction0., 10., In47:= f=x/. q1; g=y/. q 1; In48= f0.4 Out48= 0.375078 In49= g4Out49

47、= -0.0912007In50= ParametricPlotft,gt,t,0,101)0()0(22yxyxyxyxEuler方法方法lEuler方法算法方法算法l l輸入函數(shù)f(x,y)、初值y0、自變量區(qū)間端點a,b步長hl計算節(jié)點數(shù)n和節(jié)點x k l用Euler公式y(tǒng) n+1= yn+h f(xn,yn) 求數(shù)值解lEuler方法程序方法程序lClearx,y,flfx_,y_= Input函數(shù)f(x,y)=ly0=Input初值y0 =la=Input左端點a=lb=Input右端點b=lh=Input步長h=ln=(b-a)/h;lFori=0,iInterpolatingFunctionrange, 的形式給出的,其中的InterpolatingFunctionrange, 是所求的插值函數(shù)表示的數(shù)值解, range就是所求數(shù)值解的自變量范圍。l用Euler方法求初值問題 的數(shù)值解。取步長h=0.1,并在一個坐標(biāo)系中畫出數(shù)值解與準(zhǔn)確解的圖形。l解：解：執(zhí)行Euler方法程序后，在輸入的窗口中按提示分別輸入y-2x/y、1、0、1、0.1，每次輸入后用鼠標(biāo)點