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

MATLAB在數(shù)值分析中的應(yīng)用何仁斌11MATLAB在數(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é)問題●方程求根問題●解線性方程組的問題●定積分問題●常微分方程初值問題…等等在科學(xué)計(jì)算中常要遇到求解各種方程對于高次代數(shù)方程，由代數(shù)根本定理知多項(xiàng)式根的個(gè)數(shù)和方程的階相同但對超越方程就復(fù)雜的多，如果有解，其解可能是一個(gè)或幾個(gè)，也可能是無窮多個(gè)。11.1.1方程求根問題例如：高次代數(shù)方程x5

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

–cosx=0看似簡單，但難求其精確解。方程求根問題由線性代數(shù)知識可知：當(dāng)線性方程組Ax=b的系數(shù)矩陣A非奇異(即detA≠0)時(shí)，方程組有唯一解，可用克萊默法那么求解但它只適合于n很小的情況，而完全不適合于高次方程組。如用克萊默法那么求解一個(gè)n階方程組，要算n+1個(gè)n階行列式的值，總共需要n!(n-1)(n+1)次乘法。當(dāng)n充分大時(shí)，計(jì)算量是相當(dāng)驚人的解線性方程組的問題如用克萊默法那么求解一個(gè)n階方程組，要算n+1個(gè)n階行列式的值，總共需要n!(n-1)(n+1)次乘法。當(dāng)n充分大時(shí)，計(jì)算量是相當(dāng)驚人的一個(gè)20階不算太大的方程組，大約要做1021次乘法，這項(xiàng)計(jì)算即使每秒1萬億次浮點(diǎn)數(shù)乘法計(jì)算的計(jì)算機(jī)去做，也要連續(xù)工作2000萬億年才能完成。當(dāng)然這是完全沒有實(shí)際意義的，故需要尋找有效算法11.1.2解線性方程組的問題由微積分知識知，定積分的計(jì)算可以使用牛頓——萊布尼茲公式：其中F(x)為被積函數(shù)f(x)的原函數(shù)。為何要進(jìn)行數(shù)值積分？定積分問題原因之一：許多形式上很簡單的函數(shù)，例如等，它們的原函數(shù)不能用初等函數(shù)表示成有限形式。定積分問題原因之二：有些被積函數(shù)的原函數(shù)過于復(fù)雜，計(jì)算不便。例如的一個(gè)原函數(shù)是定積分問題原因之三：f(x)以離散數(shù)據(jù)點(diǎn)形式給出：xix0x1…xnyi=f(xi)y0y1…yn定積分問題對一些典型的微分方程，如可別離變量方程、一階線性方程等，有可能找出它們的一般解表達(dá)式，然后用初始條件確定表達(dá)式中的任意常數(shù)，這樣即能確定解但是對于常微分方程初值問題：那么無法求出一般解常微分方程初值問題11.2多項(xiàng)式與插值來源于實(shí)際、又廣泛用于實(shí)際。多項(xiàng)式插值的主要目的是用一個(gè)多項(xiàng)式擬合離散點(diǎn)上的函數(shù)值，使得可以用該多項(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ù)根。假設(shè)多項(xiàng)式所有系數(shù)均為實(shí)數(shù)，那么全部復(fù)根都將以共軛對的形式出現(xiàn)冪系數(shù)：在MATLAB里，多項(xiàng)式用行向量表示，其元素為多項(xiàng)式的系數(shù)，并從左至右按降冪排列。

例：

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

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

r=0.2500+1.5612i0.2500-1.5612i-1.0000所有零點(diǎn)由一個(gè)列向量給出。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=23.812576.0000polyfit:給定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.20151.4385-2.74775.4370>>poly2sym(a)ans=-403/2000*x^3+2877/2000*x^2-27477/10000*x+5437/1000多項(xiàng)式為Polyfit的第三個(gè)參數(shù)是多項(xiàng)式的階數(shù)。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)例

>>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

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

：最鄰近插值；‘spline’

：三次樣條插值；‘cubic’

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

例：在1-12的11小時(shí)內(nèi)，每隔1小時(shí)測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)每隔1/10小時(shí)的溫度值。用MATLAB作分段線性插值計(jì)算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作分段線性插值計(jì)算2024/3/14程序運(yùn)行結(jié)果：用MATLAB作分段線性插值計(jì)算2024/3/14返回

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

：最鄰近插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值;缺省時(shí)：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍。用MATLAB作插值計(jì)算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)xi處的插值結(jié)果2024/3/14例:在上取11個(gè)點(diǎn)，作三次樣條插值,觀察三次樣條插值曲線與g(x)的誤差.用MATLAB作三次樣條插值計(jì)算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作三次樣條插值計(jì)算2024/3/14返回

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

xyx1x2x3x4x5y1y2y3

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

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

‘nearest’

最鄰近插值；‘linear’

雙線性插值；

‘cubic’

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

注意：x0,y0為向量，但z0是矩陣，其列數(shù)等于x0的長度，行數(shù)等于y0的長度。例某實(shí)驗(yàn)對一根長10米的鋼軌進(jìn)行熱源的溫度傳播測試。用x表示測量點(diǎn)0:2.5:10(米)，用h表示測量時(shí)間0:30:60(秒)，用T表示測試所得各點(diǎn)的溫度(℃)。試用線性插值求出在一分鐘內(nèi)每隔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作散點(diǎn)數(shù)據(jù)的插值計(jì)算

注意：x0,y0,z0均為向量，長度相等。Method可取‘nearest’,’linear’,’cubic’;‘linear’是缺省值。z=griddata(x0,y0,z0,x,y,’method’)被插值點(diǎn)插值節(jié)點(diǎn)被插值點(diǎn)的函數(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個(gè)點(diǎn)〔xi,yi)i=1,…n,尋求一個(gè)函數(shù)〔曲線〕y=f(x),使f(x)在某種準(zhǔn)那么下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。+++++++++xyy=f(x)(xi,yi)

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

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

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

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

待定常數(shù)函數(shù)M文件擬合函數(shù)y=f(a,x)的函數(shù)M—文件各數(shù)據(jù)點(diǎn)處的誤差誤差平方和a的初值數(shù)據(jù)點(diǎn)用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ù)點(diǎn)：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,表示輸出中間結(jié)果。solve('f1(x)’,'f2(x)’,…,'fn(x)

’)X=fsolve(‘fun’,X0,options)11.4解方程組3、多項(xiàng)式方程：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特點(diǎn)：可以找出全部根。特點(diǎn)：只能求出一個(gè)特解。輸出:[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個(gè)方程的系統(tǒng)，輸入[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10]=solve(……)不但費(fèi)時(shí)，而且非常笨拙。solve可以給出結(jié)構(gòu)輸出形式。S=solve('x^2*y^2-2*x-1=0','x^2-y^2-1=0')S=x:[8x1sym]y：[8x1sym]〔給出結(jié)構(gòu)解〕輸出具體的解：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;運(yùn)行程序(test4.m)y=fsolve('nxxf',[1,1,1],1)3〕運(yùn)行結(jié)果：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：求解多項(xiàng)式方程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①簡單的微分方程。②復(fù)雜、大型的微分方程。一階微分方程:獲取解析解的方法歸類:①別離變量法；如dy/dx=x*y；②齊次方程的變換法；如dy/dx=f(y/x)③線性方程的常數(shù)變易法或公式法.……解析解MATLAB軟件實(shí)現(xiàn)解析解dsolve('eqn1','eqn2',…,'c1',…,'var1',…)微分方程組初值條件變量組注意：①y'Dy，y''D2y②自變量名可以省略，默認(rèn)變量名‘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')上述兩例的計(jì)算結(jié)果怎樣？由此得出什么結(jié)論？例③非常系數(shù)的二階微分方程例③無解析表達(dá)式！x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非線性微分方程x=[sin(t)][-sin(t)]假設(shè)欲求解的某個(gè)數(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：運(yùn)用組合的4/5階龍格-庫塔算法自變量值函數(shù)值用于設(shè)定誤差限(缺省時(shí)設(shè)定相對誤差10-3,絕對誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對誤差和絕對誤差.Matlab軟件計(jì)算數(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標(biāo)準(zhǔn)形式：y’=f(x,y)1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-函數(shù)文件中的待解方程組應(yīng)以x的分量形式寫成.2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.注意:注意1：1、建立M文件函數(shù)functionxdot=fun(t,x,y)xdot=[f1(t,x(t),y(t))；f2(t,x(t),y(t))]；2、數(shù)值計(jì)算〔執(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,藍(lán)色曲線——y〔1〕；〔原方程解〕

紅色曲線——y〔2〕；計(jì)算結(jié)果例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