MATLAB在數(shù)值分析中的應用

MATLAB在數(shù)值分析中的應用何仁斌11MATLAB在數(shù)值分析中的應用11.1為什么要學習數(shù)值分析11.2多項式與插值11.3數(shù)據(jù)擬合11.4〔非〕線性方程組11.5微分方程組11.6積分11.1為什么要學習數(shù)值分析現(xiàn)實世界的問題可以歸結為各種各樣的數(shù)學問題●方程求根問題●解線性方程組的問題●定積分問題●常微分方程初值問題…等等在科學計算中常要遇到求解各種方程對于高次代數(shù)方程，由代數(shù)根本定理知多項式根的個數(shù)和方程的階相同但對超越方程就復雜的多，如果有解，其解可能是一個或幾個，也可能是無窮多個。11.1.1方程求根問題例如：高次代數(shù)方程x5

–3x＋1=0超越方程e-x

–cosx=0看似簡單，但難求其精確解。方程求根問題由線性代數(shù)知識可知：當線性方程組Ax=b的系數(shù)矩陣A非奇異(即detA≠0)時，方程組有唯一解，可用克萊默法那么求解但它只適合于n很小的情況，而完全不適合于高次方程組。如用克萊默法那么求解一個n階方程組，要算n+1個n階行列式的值，總共需要n!(n-1)(n+1)次乘法。當n充分大時，計算量是相當驚人的解線性方程組的問題如用克萊默法那么求解一個n階方程組，要算n+1個n階行列式的值，總共需要n!(n-1)(n+1)次乘法。當n充分大時，計算量是相當驚人的一個20階不算太大的方程組，大約要做1021次乘法，這項計算即使每秒1萬億次浮點數(shù)乘法計算的計算機去做，也要連續(xù)工作2000萬億年才能完成。當然這是完全沒有實際意義的，故需要尋找有效算法11.1.2解線性方程組的問題由微積分知識知，定積分的計算可以使用牛頓——萊布尼茲公式：其中F(x)為被積函數(shù)f(x)的原函數(shù)。為何要進行數(shù)值積分？定積分問題原因之一：許多形式上很簡單的函數(shù)，例如等，它們的原函數(shù)不能用初等函數(shù)表示成有限形式。定積分問題原因之二：有些被積函數(shù)的原函數(shù)過于復雜，計算不便。例如的一個原函數(shù)是定積分問題原因之三：f(x)以離散數(shù)據(jù)點形式給出：xix0x1…xnyi=f(xi)y0y1…yn定積分問題對一些典型的微分方程，如可別離變量方程、一階線性方程等，有可能找出它們的一般解表達式，然后用初始條件確定表達式中的任意常數(shù)，這樣即能確定解但是對于常微分方程初值問題：那么無法求出一般解常微分方程初值問題11.2多項式與插值來源于實際、又廣泛用于實際。多項式插值的主要目的是用一個多項式擬合離散點上的函數(shù)值，使得可以用該多項式估計數(shù)據(jù)點之間的函數(shù)值?？蓪С鰯?shù)值積分方法，有限差分近似關注插值多項式的表達式、精度、選點效果。11.2.1關于多項式MATLAB命令一個多項式的冪級數(shù)形式可表示為：也可表為嵌套形式或因子形式

N階多項式n個根，其中包含重根和復根。假設多項式所有系數(shù)均為實數(shù)，那么全部復根都將以共軛對的形式出現(xiàn)冪系數(shù)：在MATLAB里，多項式用行向量表示，其元素為多項式的系數(shù)，并從左至右按降冪排列。

例：

被表示為>>p=[2145]>>poly2sym(p)ans=2*x^3+x^2+4*x+5Roots:多項式的零點可用命令roots求的。

例：>>r=roots(p)得到

r=0.2500+1.5612i0.2500-1.5612i-1.0000所有零點由一個列向量給出。polyval:可用命令polyval計算多項式的值。例：計算y(2.5)

>>c=[3,-7,2,1,1]；xi=2.5;yi=polyval(c,xi)yi=23.8125如果xi是含有多個橫坐標值的數(shù)組，那么yi也為與xi長度相同的向量。>>c=[3,-7,2,1,1];xi=[2.5,3];>>yi=polyval(c,xi)yi=23.812576.0000polyfit:給定n+1個點將可以唯一確定一個n階多項式。利用命令polyfit可容易確定多項式的系數(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.20151.4385-2.74775.4370>>poly2sym(a)ans=-403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000多項式為Polyfit的第三個參數(shù)是多項式的階數(shù)。m階多項式與n階多項式的乘積是d=m+n階的多項式:計算系數(shù)的MATLAB命令是:c=conv(a,b)多項式除多項式的除法滿足:其中是商,是除法的余數(shù)。多項式和可由命令deconv算出。例：[q,r]=deconv(a,b)例

>>a=[2,-5,6,-1,9];b=[3,-90,-18];>>c=conv(a,b)c=6-195432-4539-792-162>>[q,r]=deconv(c,b)q=2-56-19r=0000000>>poly2sym(c)ans=6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162

分段插值：通過插值點用折線或低次曲線連接起來逼近原曲線。MATLAB實現(xiàn)可調用內部函數(shù)。命令1interp1功能：一維數(shù)據(jù)插值〔表格查找〕。該命令對數(shù)據(jù)點之間計算內插值。它找出一元函數(shù)f(x)在中間點的數(shù)值。其中函數(shù)f(x)由所給數(shù)據(jù)決定。格式1yi=interp1(x,Y,xi)%返回插值向量yi，每一元素對應于參量xi，同時由向量x與Y的內插值決定。參量x指定數(shù)據(jù)Y的點。假設Y為一矩陣，那么按Y的每列計算。11.2.2一維插值2024/3/14Method可取：‘nearest’

：最鄰近插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值;缺省時：分段線性插值。注意：所有的插值方法都要求x是單調的，并且xi不能夠超過x的范圍。1一維插值一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結果2024/3/14

例：在1-12的11小時內，每隔1小時測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計每隔1/10小時的溫度值。用MATLAB作分段線性插值計算2024/3/14hours=1:12;temps=[589152529313022252724];h=1:.1:12;t=interp1(hours,temps,h)plot(hours,temps,'+',h,t)title('線性插值下的溫度曲線'),xlabel('Hour'),ylabel('DegreesCelsius')用MATLAB作分段線性插值計算2024/3/14程序運行結果：用MATLAB作分段線性插值計算2024/3/14返回

2024/3/14Method可取：‘nearest’

：最鄰近插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值;缺省時：分段線性插值。注意：所有的插值方法都要求x是單調的，并且xi不能夠超過x的范圍。用MATLAB作插值計算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結果2024/3/14例:在上取11個點，作三次樣條插值,觀察三次樣條插值曲線與g(x)的誤差.用MATLAB作三次樣條插值計算2024/3/14x0=linspace(-5,5,11);y0=1./(1+x0.^2);x=linspace(-5,5,100);y=interp1(x0,y0,x,'spline');x1=linspace(-5,5,100);y1=1./(1+x1.^2);plot(x1,y1,'k',x0,y0,'+',x,y,'r');M文件yangcha.m用MATLAB作三次樣條插值計算2024/3/14返回

2024/3/14二維插值的提法mn個節(jié)點(xi,xj,zij)(i=1,2,…,m;j=1,2,…,n)其中xi,yj互不相同，不妨設a=x1<x2<…<xm=bc=y1<y2<…<yn=d求任一插值點(x*,y*)(≠(xi,yj))處的插值Z*第一種〔網(wǎng)格節(jié)點〕11.2.3二維插值2024/3/14二維插值的提法

xyx1x2x3x4x5y1y2y3

(x*,y*)2024/3/14第二種〔散亂節(jié)點〕n個節(jié)點其中互不相同，求任一插值點處的插值二維插值的提法2024/3/14二維插值方法的根本思路

構造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即或返回2024/3/14Method可?。?/p>

‘nearest’

最鄰近插值；‘linear’

雙線性插值；

‘cubic’

雙三次插值；缺省時,雙線性插值。二維插值:已有程序z=interp2(x0,y0,z0,x,y,’method’)被插值點插值方法插值節(jié)點被插值點的函數(shù)值網(wǎng)格節(jié)點數(shù)據(jù)的插值

注意：x0,y0為向量，但z0是矩陣，其列數(shù)等于x0的長度，行數(shù)等于y0的長度。例某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度(℃)。試用線性插值求出在一分鐘內每隔20秒、鋼軌每隔1米處的溫度TI。命令如下：x=0:2.5:10;h=0:30:60;T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)2024/3/14用MATLAB作散點數(shù)據(jù)的插值計算

注意：x0,y0,z0均為向量，長度相等。Method可取‘nearest’,’linear’,’cubic’;‘linear’是缺省值。z=griddata(x0,y0,z0,x,y,’method’)被插值點插值節(jié)點被插值點的函數(shù)值x=rand(100,1)*4-2;y=rand(100,1)*4-2;z=x.*exp(-x.^2-y.^2);ti=-2:.25:2;[XI,YI]=meshgrid(ti,ti);ZI=griddata(x,y,z,XI,YI);mesh(XI,YI,ZI),holdplot3(x,y,z,'o'),holdoff11.3曲線擬合一組〔二維〕數(shù)據(jù)，即平面上n個點〔xi,yi)i=1,…n,尋求一個函數(shù)〔曲線〕y=f(x),使f(x)在某種準那么下與所有數(shù)據(jù)點最為接近，即曲線擬合得最好。+++++++++xyy=f(x)(xi,yi)

ii為點〔xi,yi)與曲線y=f(x)的距離問題的數(shù)學模型步驟：1)選定一類函數(shù)f(x，a1，a2，…，am)〔1〕其中a1,a2,…am為待定常數(shù)。2)確定參數(shù)a1,a2,…am準那么(最小二乘準那么〕：使n個點〔xi,yi)與曲線y=f(x，a1，a2，…，am)的距離i的平方和最小。曲線擬合的根本原理記

問題歸結為:求a1,a2,…am使J(a1,a2,…am)最小。這樣的擬合稱為最小二乘擬合。問題的數(shù)學模型曲線擬合的根本原理最小二乘擬合函數(shù)的選取+++++++++++++++f=a1+a2x將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過直觀判斷確定f：f=a1+a2x+a3x2f=a1+a2x+a3x22.通過機理分析建立數(shù)學模型來確定f。+++++++++++++++f=a1+a2/xf=a1exp(a2x)f=a1exp(a2x)最小二乘擬合函數(shù)的選取多項式擬合：作多項式f(x)=a1xm+…+amx+am+1函數(shù)擬合,可利用已有程序polyfit,其調用格式為:a=polyfit(x,y,m)用MATLAB作最小二乘擬合數(shù)據(jù)點系數(shù)擬合多項式次數(shù)例.熱敏電阻數(shù)據(jù)：溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032求電阻R隨溫度t的變化規(guī)律。1.多項式擬合用MATLAB作最小二乘擬合2)用命令polyfit(x,y,m)作最小二乘擬合3)編寫MATLAB程序dianzu.m，并運行得到：a1=3.3940,a2=702.49181)選取擬合函數(shù)R=a1t+a21.多項式擬合用MATLAB作最小二乘擬合t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')M文件dianzu.m1.多項式擬合用MATLAB作最小二乘擬合2.曲線擬合：

作一般的最小二乘曲線擬合，可利用已有程序lsqcurvefit,其調用格式為：

[a,resnorm,residual]=lsqcurvefit(‘f’,a0,x,y)

待定常數(shù)函數(shù)M文件擬合函數(shù)y=f(a,x)的函數(shù)M—文件各數(shù)據(jù)點處的誤差誤差平方和a的初值數(shù)據(jù)點用MATLAB作最小二乘擬合xdata=[0:0.1:2];ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]例：用函數(shù)y(x)=a1*exp(-a2*x)+a3*exp(-a4*x)擬合以下數(shù)據(jù)點：1.用命令lsqcurvefit(fun,a0,x,y)2.曲線擬合用MATLAB作最小二乘擬合f=@(a,x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x);a0=[1,1,1,0];[a,resnorm,residual,flag,output]=lsqcurvefit(f,a0,xdata,ydata)xi=linspace(0,2,200);yi=f(a,xi);plot(xdata,ydata,'ro',xi,yi)xlabel('x'),ylabel('y=f(x)'),title('nonlinearcurvefitting')2.曲線擬合用MATLAB作最小二乘擬合1、方程(組)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)求解方程(組)的MATLAB語句2、方程(組)，f1(x)=0,…,fn(x)=0,x=(x1,…,xn)fun.mfunctionf=fun(x)f(1)=f1(x);……f(n)=fn(x)初值1）可以省略。2）options=1,表示輸出中間結果。solve('f1(x)’,'f2(x)’,…,'fn(x)

’)X=fsolve(‘fun’,X0,options)11.4解方程組3、多項式方程：amxm+am-1xm-1+…+a0=0p=[am,am-1,…,a0]；roots(p)4、線性方程組：AX=b

其中A是m×n階矩陣，b是m維向量。x=A\borx=inv(A)*b特點：可以找出全部根。特點：只能求出一個特解。輸出:[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]①單變量方程solve()語句的用法例1：求解方程ax2+bx+c=0輸入: x=solve('a*x^2+b*x+c')或solve('a*x^2+b*x+c=0')1〕符號解例2：解方程：x3-2x2=x-1解：

s=solve('x^3-2*x^2=x-1')double(s)2〕數(shù)值解solve()語句的用法s=7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))+((108^(1/2)*7*i)/108+7/54)^(1/3)+2/32/3-((108^(1/2)*7*i)/108+7/54)^(1/3)/2+(3^(1/2)*(7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))-((108^(1/2)*7*i)/108+7/54)^(1/3))*i)/2-7/(18*(108^(1/2)*((7*i)/108)+7/54)^(1/3))2/3-((108^(1/2)*7*i)/108+7/54)^(1/3)/2-(3^(1/2)*(7/(9*(108^(1/2)*((7*i)/108)+7/54)^(1/3))-((108^(1/2)*7*i)/108+7/54)^(1/3))*i)/2-7/(18*(108^(1/2)*((7*i)/108)+7/54)^(1/3))ans=2.2470+0.0000i0.5550-0.0000i-0.8019-0.0000i輸入：[x,y]=solve('x^2*y^2,x-(y/2)-b')輸出：x=[0]，y=[-2*b][0]，[-2*b]〔符號解〕[b]，[0][b]，[0]v=[x,y]②多變量方程組例4solve()語句的用法如果有10個方程的系統(tǒng)，輸入[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]=solve(……)不但費時，而且非常笨拙。solve可以給出結構輸出形式。S=solve('x^2*y^2-2*x-1=0','x^2-y^2-1=0')S=x:[8x1sym]y：[8x1sym]〔給出結構解〕輸出具體的解：S.x,S.y,s1=[S.x(2),S.y(2)]；〔取出解空間的第二組解〕M=[S.x,S.y]；〔創(chuàng)立解矩陣〕例5solve()語句的用法例6：求解方程組fsolve()語句的用法解①輸入：[x,y,z]=solve('sin(x)+y^2+log(z)-7=0','3*x+2^y-z^3+1=0','x+y+z-5=0','x','y','z')x=0.59905375664056731520568183824539y=2.3959314023778168490940003756591z=2.0050148409816158357003177860955輸出：解②：1〕建立方程組的M-函數(shù)文件(nxxf.m)functioneq=nxxf(x)eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;eq(2)=3*x(1)+2^x(2)-x(3)^3+1;eq(3)=x(1)+x(2)+x(3)-5;運行程序(test4.m)y=fsolve('nxxf',[1,1,1],1)3〕運行結果：OptimizationTerminatedSuccessfullyy=0.59902.39592.0050fsolve()語句的用法輸出:

-1.2131-0.9017+0.5753i-0.9017-0.5753i-0.2694+0.9406i-0.2694-0.9406i0.4168+0.8419i0.4168-0.8419i0.8608+0.3344i0.8608-0.3344i例7：求解多項式方程x9+x8+1=0roots()語句的用法輸入:p=[1,1,0,0,0,0,0,0,0,1];roots(p)例8:AX=b,A\b和inv(A)*b語句的用法解:輸入: A=[123;456;789];b=[6;14;-3];x1=A\b,x2=inv(A)*b輸出: 警告:系統(tǒng)的秩缺乏.解不唯一.特殊情形：1、微分方程的一般形式：F(x,y,y’,…,y(n))=0

隱式或y(n)=f(x,y,y’,…,y(n-1))

顯式2、一階微分方程組的一般形式：n階1階初始條件：y(x0)=y011.5解微分方程微分方程解的形式①解析解y=f(x)②數(shù)值解(xi,yi)③圖形解xyo①簡單的微分方程。②復雜、大型的微分方程。一階微分方程:獲取解析解的方法歸類:①別離變量法；如dy/dx=x*y；②齊次方程的變換法；如dy/dx=f(y/x)③線性方程的常數(shù)變易法或公式法.……解析解MATLAB軟件實現(xiàn)解析解dsolve('eqn1','eqn2',…,'c1',…,'var1',…)微分方程組初值條件變量組注意：①y'Dy，y''D2y②自變量名可以省略，默認變量名‘t’。例①輸入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')輸出：y=tan(t-C1)〔通解，一簇曲線〕y1=tan(x+1/4*pi)〔特解，一條曲線〕例②常系數(shù)的二階微分方程y=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')輸入：x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')上述兩例的計算結果怎樣？由此得出什么結論？例③非常系數(shù)的二階微分方程例③無解析表達式！x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非線性微分方程x=[sin(t)][-sin(t)]假設欲求解的某個數(shù)值解，如何求解？t=pi/2;eval(x)輸入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例④輸出：〔li3.m〕數(shù)值解1、歐拉法2、龍格—庫塔法數(shù)值求解思想：(變量離散化)

研究常微分方程的數(shù)值解法是十分必要的。[t，x]=solver〔’f’,ts,x0,options〕ode23

ode45ode113ode15sode23s由待解方程寫成的m-函數(shù)文件ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫塔算法ode45：運用組合的4/5階龍格-庫塔算法自變量值函數(shù)值用于設定誤差限(缺省時設定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設定的相對誤差和絕對誤差.Matlab軟件計算數(shù)值解f=@(x,y)-y+x+1;[x,y]=ode23(f,[0,1],1)plot(x,y,'r');holdonezplot('x+exp(-x)',[0,1])%精確解例1

y’=-y+x+1，y(0)=1標準形式：y’=f(x,y)1、在解n個未知函數(shù)的方程組時，x0和x均為n維向量，m-函數(shù)文件中的待解方程組應以x的分量形式寫成.2、使用Matlab軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組.注意:注意1：1、建立M文件函數(shù)functionxdot=fun(t,x,y)xdot=[f1(t,x(t),y(t))；f2(t,x(t),y(t))]；2、數(shù)值計算〔執(zhí)行以下命令〕[t,x,y]=ode23(‘fun',[t0,tf],[x0,y0])注意：執(zhí)行命令不能寫在M函數(shù)文件中。xd(1)=f1(t,x(t),y(t));xd(2)=f2(t,x(t),y(t));xdot=xd’;%列向量例如：令注意2：functionxdot=fun1(t,x,y)〔fun1.m)xdot=[f(t,x(t),y(t))；x(t)]；[t,x,y]=ode23(‘fun1',[t0,tf],[x0,y0])M-文件函數(shù)如何寫呢？注意：y(t)是原方程的解。x(t)只是中間變量。如果方程形式是：z’’’=f(t,z,z’’)?例2Vanderpol方程：令y1=x(t),y2=x’(t);

該方程是否有解析解？編寫程序如下:

f=@(t,y)[y(2);(1-y(1)^2)*y(2)-y(1)];[t,y]=ode23(f,[0,20],[3,0]);y1=y(:,1);%原方程的解y2=y(:,2);plot(t,y1,t,y2,'--')%y1(t),y2(t)曲線圖

pause,plot(y1,y2),%相軌跡圖，即y2(y1)曲線grid,藍色曲線——y〔1〕；〔原方程解〕

紅色曲線——y〔2〕；計算結果例3考慮Lorenz模型：其中參數(shù)β=8/3，σ=10，ρ=28解：1〕編寫匿名函數(shù)；2)數(shù)值求解并畫三維空間的相平面軌線。f=@(t,x)[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x;x0=[000.1]';[t,x]=ode45(f,[0,10],x0);pl