1、MATLAB在數(shù)值分析中的應(yīng)用 何仁斌 11 MATLAB在數(shù)值分析中的應(yīng)用11.1為什么要學(xué)習(xí)數(shù)值分析11.2 多項(xiàng)式與插值11.3 數(shù)據(jù)擬合11.4 （非）線性方程組11.5 微分方程組11.6 積分11.1為什么要學(xué)習(xí)數(shù)值分析現(xiàn)實(shí)世界的問題可以歸結(jié)為各種各樣的數(shù)學(xué)問題現(xiàn)實(shí)世界的問題可以歸結(jié)為各種各樣的數(shù)學(xué)問題 方程求根問題方程求根問題 解線性方程組的問題解線性方程組的問題 定積分問題定積分問題 常微分方程初值問題常微分方程初值問題 等等等等在科學(xué)計(jì)算中常要遇到求解各種方程在科學(xué)計(jì)算中常要遇到求解各種方程對于高次代數(shù)方程，由代數(shù)基本定理知多對于高次代數(shù)方程，由代數(shù)基本定理知多項(xiàng)式根的個(gè)數(shù)和

2、方程的階相同項(xiàng)式根的個(gè)數(shù)和方程的階相同但對超越方程就復(fù)雜的多，如果有解，其但對超越方程就復(fù)雜的多，如果有解，其解可能是一個(gè)或幾個(gè)，也可能是無窮多解可能是一個(gè)或幾個(gè)，也可能是無窮多個(gè)。個(gè)。11.1.1 方程求根問題例如：例如：高次代數(shù)方程高次代數(shù)方程 x5 3x1 = 0超越方程超越方程 e-x cosx = 0看似簡單，但難求其精確解。看似簡單，但難求其精確解。11.1.1方程求根問題由線性代數(shù)知識可知：當(dāng)線性方程組由線性代數(shù)知識可知：當(dāng)線性方程組Ax = b的系數(shù)矩陣的系數(shù)矩陣A非奇異非奇異(即即detA0)時(shí)，方程時(shí)，方程組有唯一解，可用克萊默法則求解組有唯一解，可用克萊默法則求解但它只適

3、合于但它只適合于n很小的情況，而完全不適很小的情況，而完全不適合于高次方程組。合于高次方程組。如用克萊默法則求解一個(gè)如用克萊默法則求解一個(gè)n階方程組，要階方程組，要算算n+1個(gè)個(gè)n階行列式的值，總共需要階行列式的值，總共需要n!(n-1)(n+1)次乘法。當(dāng)次乘法。當(dāng)n充分大時(shí)，計(jì)算量是充分大時(shí)，計(jì)算量是相當(dāng)驚人的相當(dāng)驚人的11.1.2解線性方程組的問題如用克萊默法則求解一個(gè)如用克萊默法則求解一個(gè)n階方程組，要階方程組，要算算n+1個(gè)個(gè)n階行列式的值，總共需要階行列式的值，總共需要n!(n-1)(n+1)次乘法。當(dāng)次乘法。當(dāng)n充分大時(shí)，計(jì)算量是充分大時(shí)，計(jì)算量是相當(dāng)驚人的相當(dāng)驚人的一個(gè)一個(gè)20

4、階不算太大的方程組，大約要做階不算太大的方程組，大約要做1021次乘法，這項(xiàng)計(jì)算即使每秒次乘法，這項(xiàng)計(jì)算即使每秒1萬億次萬億次浮點(diǎn)數(shù)乘法計(jì)算的計(jì)算機(jī)去做，也要連續(xù)浮點(diǎn)數(shù)乘法計(jì)算的計(jì)算機(jī)去做，也要連續(xù)工作工作2000萬億年才能完成。萬億年才能完成。當(dāng)然這是完全沒有實(shí)際意義的，故需要尋當(dāng)然這是完全沒有實(shí)際意義的，故需要尋找有效算法找有效算法11.1.2 解線性方程組的問題由微積分知識知，定積分的計(jì)算可以使用由微積分知識知，定積分的計(jì)算可以使用牛頓牛頓萊布尼茲公式：萊布尼茲公式：其中其中F(x)為被積函數(shù)為被積函數(shù)f(x)的原函數(shù)。的原函數(shù)。為何要進(jìn)行數(shù)值積分？為何要進(jìn)行數(shù)值積分？( )( )( )

5、baf x dxF bF a11.1.3定積分問題原因之一：許多形式上很簡單的函數(shù)，例原因之一：許多形式上很簡單的函數(shù)，例如如等，它們的原函數(shù)不能用初等函數(shù)表示成等，它們的原函數(shù)不能用初等函數(shù)表示成有限形式。有限形式。2sinsin,bbaaxx dxdxx11.1.3定積分問題原因之二：有些被積函數(shù)的原函數(shù)過于復(fù)原因之二：有些被積函數(shù)的原函數(shù)過于復(fù)雜，計(jì)算不便。例如雜，計(jì)算不便。例如的一個(gè)原函數(shù)是的一個(gè)原函數(shù)是22( )23f xxx3222233239( )ln( 223)41616 2xxxxF xxx11.1.3定積分問題原因之三：原因之三：f(x)以離散數(shù)據(jù)點(diǎn)形式給出：以離散數(shù)據(jù)點(diǎn)形

6、式給出：xix0 x1xnyi = f(xi)y0y1yn11.1.3定積分問題對一些典型的微分方程，如可分離變量方對一些典型的微分方程，如可分離變量方程、一階線性方程等，有可能找出它們的程、一階線性方程等，有可能找出它們的一般解表達(dá)式，然后用初始條件確定表達(dá)一般解表達(dá)式，然后用初始條件確定表達(dá)式中的任意常數(shù)，這樣即能確定解式中的任意常數(shù)，這樣即能確定解但是對于常微分方程初值問題：但是對于常微分方程初值問題：則無法求出一般解則無法求出一般解2(0)0yyxy11.1.4常微分方程初值問題11.2 多項(xiàng)式與插值 來源于實(shí)際、又廣泛用于實(shí)際。 多項(xiàng)式插值的主要目的是用一個(gè)多項(xiàng)式擬合離散點(diǎn)上的函數(shù)值

7、，使得可以用該多項(xiàng)式估計(jì)數(shù)據(jù)點(diǎn)之間的函數(shù)值。 可導(dǎo)出數(shù)值積分方法，有限差分近似 關(guān)注插值多項(xiàng)式的表達(dá)式、精度、選點(diǎn)效果。11.2.1 關(guān)于多項(xiàng)式MATLAB命令 一個(gè)多項(xiàng)式的冪級數(shù)形式可表示為： 也可表為嵌套形式 或因子形式 N階多項(xiàng)式n個(gè)根，其中包含重根和復(fù)根。若多項(xiàng)式所有系數(shù)均為實(shí)數(shù)，則全部復(fù)根都將以共軛對的形式出現(xiàn) 1121nnnnyc xc xc xc1231()nnyc xc xc xc xc112()()()nyc xrxrxr 冪系數(shù)：在MATLAB里，多項(xiàng)式用行向量表示，其元素為多項(xiàng)式的系數(shù)，并從左至右按降冪排列。 例： 被表示為 p=2 1 4 5 poly2sym(p) a

8、ns = 2*x3+x2+4*x+5 Roots: 多項(xiàng)式的零點(diǎn)可用命令roots求的。 例： r=roots(p) 得到 r = 0.2500 + 1.5612i 0.2500 - 1.5612i -1.0000 所有零點(diǎn)由一個(gè)列向量給出。32245yxxx polyval: 可用命令polyval計(jì)算多項(xiàng)式的值。 例： 計(jì)算y(2.5) c=3,-7,2,1,1； xi=2.5; yi=polyval(c,xi) yi = 23.8125如果xi是含有多個(gè)橫坐標(biāo)值的數(shù)組，則yi也為與xi長度相同的向量。 c=3,-7,2,1,1; xi=2.5,3; yi=polyval(c,xi)yi

9、= 23.8125 76.00004323721yxxxx polyfit:給定n+1個(gè)點(diǎn)將可以唯一確定一個(gè)n階多項(xiàng)式。利用命令polyfit可容易確定多項(xiàng)式的系數(shù)。 例： x=1.1,2.3,3.9,5.1; y=3.887,4.276,4.651,2.117; a=polyfit(x,y,length(x)-1)a = -0.2015 1.4385 -2.7477 5.4370 poly2sym(a) ans = -403/2000*x3+2877/2000*x2-27477/10000*x+5437/1000 多項(xiàng)式為Polyfit的第三個(gè)參數(shù)是多項(xiàng)式的階數(shù)。320.20151.4385

10、2.74775.4370yxxx m階多項(xiàng)式與n階多項(xiàng)式的乘積是d=m+n階的多項(xiàng)式:計(jì)算 系數(shù)的MATLAB命令是:c=conv(a,b) 多項(xiàng)式 除多項(xiàng)式 的除法滿足:其中 是商, 是除法的余數(shù)。多項(xiàng)式 和 可由命令deconv算出。 例：q, r=deconv(a,b)1121ddcabddyy yc xc xc xccybyayaqbryy yyqyryqyry1121mmammya xa xa xa1121nnbnnyb xb xb xb 例 a=2,-5,6,-1,9; b=3,-90,-18; c=conv(a,b)c = 6 -195 432 -453 9 -792 -162

11、q,r=deconv(c,b)q = 2 -5 6 -1 9r = 0 0 0 0 0 0 0 poly2sym(c) ans = 6*x6-195*x5+432*x4-453*x3+9*x2-792*x-162 分段插值：通過插值點(diǎn)用折線或低次曲線連接起來逼近原曲線。 MATLAB實(shí)現(xiàn) 可調(diào)用內(nèi)部函數(shù)。 命令1 interp1 功能 ： 一維數(shù)據(jù)插值（表格查找）。該命令對數(shù)據(jù)點(diǎn)之間計(jì)算內(nèi)插值。它找出一元函數(shù)f(x)在中間點(diǎn)的數(shù)值。其中函數(shù)f(x)由所給數(shù)據(jù)決定。 格式1 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素對應(yīng)于參量xi，同時(shí)由向量x與Y的內(nèi)插值決定。參量x

12、指定數(shù)據(jù)Y的點(diǎn)。若Y為一矩陣，則按Y的每列計(jì)算。11.2.2 一維插值2022-5-13Method可?。嚎扇。簄earest ：最鄰近插值；：最鄰近插值；spline ：三次樣條插值；：三次樣條插值；cubic ： 立方插值立方插值;缺省時(shí)：缺省時(shí)： 分段線性插值。分段線性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x是單調(diào)的，并且是單調(diào)的，并且xi不能夠超過不能夠超過x的范圍。的范圍。1 1 一維插值一維插值一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xi處的插處的插值結(jié)果值結(jié)果2022-5-1

13、32022-5-132022-5-13程序運(yùn)行結(jié)果：2022-5-132022-5-13Method可?。嚎扇。簄earest ：最鄰近插值；：最鄰近插值；spline ：三次樣條插值；：三次樣條插值；cubic ： 立方插值立方插值;缺省時(shí)：缺省時(shí)： 分段線性插值。分段線性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x是單調(diào)的，并且是單調(diào)的，并且xi不能夠超過不能夠超過x的范圍。的范圍。用用MATLABMATLAB作插值計(jì)算作插值計(jì)算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xi處的插處的插值結(jié)

14、果值結(jié)果2022-5-13例例: :在在21( ),551g xxx上取上取1111個(gè)點(diǎn)，作三次樣條插值個(gè)點(diǎn)，作三次樣條插值, , 觀觀察三次樣條插值曲線與察三次樣條插值曲線與g(x)g(x)的誤差的誤差. .用用MATLABMATLAB作三次樣條插值計(jì)算作三次樣條插值計(jì)算2022-5-13M文件文件yangcha.m用用MATLABMATLAB作三次樣條插值計(jì)算作三次樣條插值計(jì)算2022-5-132022-5-13二二維維插插值值的的提提法法已知已知 m n個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn)(xi, xj, zij) ( i=1, 2, ,m; j=1, 2, , n )其中其中xi, yj互不相同，不妨設(shè)互不相

15、同，不妨設(shè)a=x1x2xm=bc=y1y2 y1=dblquad(x,y) sqrt(104-x.2).*(x.2+y.2 y2=quad2d(x,y) sqrt(104-x.2),-100,100,.(x) -sqrt(104-x.2),(x) sqrt(104-x.2)y2 = 2.6667e+006后者運(yùn)行時(shí)間更少一般區(qū)域二重積分計(jì)算下列二重積分計(jì)算下列二重積分 可用符號積分計(jì)算: syms x y int(int(x*y,y,sin(x),cos(x),1,2)ans =cos(4)/8 - cos(2)/8 - sin(2)/4 + sin(4)/2 vpa(ans,10)ans =-0.63541270242cos( )1sin( )xxxydydx 一般區(qū)域二重積分計(jì)算下列二重積分計(jì)算下列二重積分 quad2d(x,y) x.*y,1,2,(x) sin(x),(x) cos(x)ans = -0.63542cos( )1sin( )xxxydydx 一般區(qū)域二重積分計(jì)算下列二重積分計(jì)算下列二重積分function f1=myfun1(x) f1=zeros(