MATLAB編程基礎(chǔ)第6講-插值、擬合與初值常微分方程的求解

MATLAB編程根底

之插值、擬合與初值常微分方程的求解第六講梁丙臣5/30/20241多項(xiàng)式的擬合需求：實(shí)驗(yàn)數(shù)據(jù)總結(jié)為規(guī)律曲線等等最小二乘法：所用曲線限定為多項(xiàng)式，在數(shù)據(jù)點(diǎn)上擬合值和函數(shù)值有最小誤差平方和。P=polyfit(x,y,n)：把自變量x和函數(shù)值y擬合成n階多項(xiàng)式，向量為p。5/30/20242例3-35x=(0:0.1:1)';y=[-0.41.93.26.27.17.37.79.69.59.312]';%給出一組11個(gè)點(diǎn)數(shù)據(jù)y2=polyfit(x,y,2)%計(jì)算2階擬合的多項(xiàng)式向量

x1=0:0.01:1;f2=polyval(y2,x1);%2階擬合曲線在各點(diǎn)的函數(shù)值y10=polyfit(x,y,10)%計(jì)算10階擬合的多項(xiàng)式向量

f10=polyval(y10,x1);%10階擬合曲線在各點(diǎn)的函數(shù)值plot(x,y,'o',x1,f2,':',x1,f10,'k')5/30/20243例3-36x=(0:0.1:2.5)';%給出一組數(shù)據(jù)，為誤差函數(shù)的一個(gè)區(qū)間y=erf(x);p=polyfit(x,y,6)%計(jì)算該區(qū)間內(nèi)6階擬合多項(xiàng)式的向量

x1=(0:0.1:5)';%將區(qū)間增長一倍y1=erf(x1);%計(jì)算誤差函數(shù)在新區(qū)間內(nèi)的函數(shù)值f=polyval(p,x1);%計(jì)算6階擬合曲線在新區(qū)間內(nèi)的取值plot(x1,y1,'o',x1,f,'-')%繪圖，比較前后區(qū)間內(nèi)的曲線擬合效果，如圖3-7axis([0502])5/30/20244例：擬合以下數(shù)據(jù)x0.51.01.52.02.53.0y1.752.453.814.808.008.60x=[0.51.01.52.02.53.0];y=[1.752.453.814.808.008.60];a=polyfit(x,y,2)a=0.49001.25010.8560x1=[0.5:0.05:3.0];y1=a(3)+a(2)*x1+a(1)*x1.^2;plot(x,y,'*')holdonplot(x1,y1,'-r')5/30/20245例：根據(jù)經(jīng)驗(yàn)公式y(tǒng)=a+bx2，擬合如下數(shù)據(jù)：xi1925313844yi19.032.349.073.398.8x=[1925313844];y=[19.032.349.073.398.8];x1=x.^2;x1=[ones(5,1),x1'];ab=x1\y'ab=0.59370.0506x0=[19:0.2:44];y0=ab(1)+ab(2)*x0.^2;plot(x,y,'o',x0,y0,'-r')5/30/20246多項(xiàng)式的插值插值的定義——是對某些集合給定的數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法。當(dāng)不能很快地求出所需中間點(diǎn)的函數(shù)時(shí)，插值是一個(gè)非常有價(jià)值的工具。Matlab提供了一維、二維、三次樣條等許多插值選擇5/30/202471.一維插值yi=interp1(x,Y,xi,method)：求同維數(shù)據(jù)x和Y，運(yùn)用method指定的方法計(jì)算插值點(diǎn)xi處的數(shù)值yi。method有如下四種：nearest最近點(diǎn)差值，取與最近的數(shù)據(jù)點(diǎn)的值linear線性差值，直線連接數(shù)據(jù)點(diǎn)，插值點(diǎn)值位于直線上spline樣條插值，用三次樣條曲線通過數(shù)據(jù)點(diǎn)，根據(jù)曲線進(jìn)行差值cubic立方差值，用三次曲線擬合并通過數(shù)據(jù)點(diǎn)5/30/20248例3-37x=0:10;%給出基準(zhǔn)數(shù)據(jù)y=sin(x);xi=0:.25:10;%在兩個(gè)基準(zhǔn)數(shù)據(jù)點(diǎn)間插入3個(gè)點(diǎn)yi1=interp1(x,y,xi,'*nearest');yi2=interp1(x,y,xi,'*linear');yi3=interp1(x,y,xi,'*spline');yi4=interp1(x,y,xi,'*cubic');%分別用四種方法求中間插值plot(x,y,'o',xi,yi1,'r:',xi,yi2,'g-',xi,yi3,'k.-',xi,yi4,'m--')%繪圖，并標(biāo)注各曲線代表的插值方法，如圖3-8legend('原始數(shù)據(jù)','最近點(diǎn)插值','線性插值','樣條插值','立方插值')5/30/202492.二維插值ZI=interp2(X,Y,Z,XI,YI,method)：同維數(shù)據(jù)X,Y和Z，運(yùn)用method指定的方法計(jì)算插值點(diǎn)(XI,YI)處的數(shù)值ZI。5/30/202410例3-38axis([-33-33-520])%確定坐標(biāo)軸的范圍[X,Y]=meshgrid(-3:.5:3);%生成網(wǎng)格Z=peaks(X,Y);%計(jì)算peaks函數(shù)值figure(1)mesh(X,Y,Z)%繪制原始數(shù)據(jù)曲面圖，如圖3-9

5/30/202411[XI,YI]=meshgrid(-3:.125:3);%確定插值點(diǎn)ZI1=interp2(X,Y,Z,XI,YI,'*nearest');%計(jì)算用nearest方法所得二維插值figure(2)mesh(XI,YI,ZI1)%繪制最近點(diǎn)插值法曲面圖，如圖3-10ZI2=interp2(X,Y,Z,XI,YI,'*linear');%計(jì)算用linear法所得二維插值figure(3)mesh(XI,YI,ZI2)%繪制雙線性插值法曲面圖，如圖3-115/30/202412ZI3=interp2(X,Y,Z,XI,YI,'*spline');%計(jì)算用spline法所得二維插值figure(4)mesh(XI,YI,ZI3)%繪制雙樣條插值法曲面圖，如圖3-12ZI4=interp2(X,Y,Z,XI,YI,'*cubic');%計(jì)算用cubic法所得二維插值figure(5)mesh(XI,YI,ZI4)%繪制雙立方插值法曲面圖，如圖3-135/30/2024133.9常微分方程初值問題的數(shù)值解法龍格－庫塔法簡介龍格－庫塔法的實(shí)現(xiàn)基于龍格－庫塔法，MATLAB提供了求常微分方程數(shù)值解的函數(shù)，一般調(diào)用格式為：[t,y]=ode23('fname',tspan,y0)[t,y]=ode45('fname',tspan,y0)其中fname是定義f(t,y)的函數(shù)文件名，該函數(shù)文件必須返回一個(gè)列向量。tspan形式為[t0,tf],表示求解區(qū)間。y0是初始狀態(tài)列向量。t和y分別給出時(shí)間向量和相應(yīng)的狀態(tài)向量。5/30/202414一階常微分方程的初值問題求解例3-39，求微分方程組在區(qū)間[012]內(nèi)的數(shù)值解，且滿足初始條件5/30/202415例3-39%首先編寫方程函數(shù)，名為rigid.mfunctiondy=rigid(t,y)dy=zeros(3,1);%生成3×1的矩陣dy(1)=y(2)*y(3);%第一個(gè)元素是dy(2)=-y(1)*y(3);%第二個(gè)元素是dy(3)=-0.51*y(1)*y(2);%第三個(gè)元素是

5/30/202416例3-39〔續(xù)前〕tspan=[012];y0=[011];%在同一目錄下，計(jì)算方程數(shù)值解。輸入時(shí)間區(qū)間和初始條件[t,Y]=ode45('rigid',tspan,y0);%采用ode45算法求解方程，options為默認(rèn)值plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')%繪制計(jì)算結(jié)果并標(biāo)注，如圖3-14legend('Y1','Y2','Y3')5/30/202417高階微分方程的初值問題求解高階微分方程：把高階微分方程轉(zhuǎn)換為一階方程組，初始條件也做相應(yīng)的替換。通常令：上式改寫為：5/30/202418例3-40，求方程初始條件在時(shí)間區(qū)間[012]內(nèi)的數(shù)值解把高階微分轉(zhuǎn)化低階微分形式，令：原方程變?yōu)椋撼跏紬l件：5/30/202419例3-40訓(xùn)練任務(wù)：請編制名為odet3.m方程函數(shù)tspan=[012];%在同一目錄下，計(jì)算方程數(shù)值解。給出時(shí)間區(qū)間和初始值y0=[-110];[t,Y]=ode45('odet3',tspan,y0);%采用ode45算法plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')%繪制計(jì)算結(jié)果并標(biāo)注，如圖3-15legend('Y1','Y2','Y3')5/30/202420%首先編寫方程函數(shù)，名為odet3.mfunctiondy=odet3(t,y)dy=zeros(3,1);dy(1)=y(2);dy(2)=y(3);dy(3)=y(1).^(1/2)-2*y(3).^2*(y(2)-4);5/30/202421例設(shè)有初值問題，試求其數(shù)值解，并與精確解相比較(精確解為y(t)=)。(1)建立函數(shù)文件funt.m。functionyp=funt(t,y)yp=(y^2-t-2)/4/(t+1);(2)求解微分方程。t0=0;tf=10;y0=2;[t,y]=ode23('funt',[t0,tf],y0);%求數(shù)值解y1=sqrt(t+1)+1;%求精確解t'y'y1'y為數(shù)值解，y1為精確值，顯然兩者近似。5/30/202422例：x+(x2-1)x+x=0為方便令x1=x，x2=x分別對x1,x2求一階導(dǎo)數(shù)，整理后寫成一階微分方程組形式

x1=x2x2=x2(1-x12)-x1建立m文件解微分方程······5/30/202423建立m文件functionxdot=wf(t,x)xdot=zeros(2,1)xdot(1)=x(2)xdot(2)=x(2)*(1-x(1)^2)-x(1)給定區(qū)間、初始值;求解微分方程t0=0;tf=20;x0=[00.25]';[t,x]=ode23('wf',t0,tf,x0)figure(1),plot(t,x)；figure(2),plot(x(:,1),x(:,2))5/30/2024245/30/202425訓(xùn)練1a=[34-7-12;5-742;108-5;-65-210];b=[4-39-8]';c=a\b5/30/202426訓(xùn)練2訓(xùn)練3算出x3+2x2+x+1=0的根functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);[S,n]=quad('fesin',0,3*pi)S=0.9008n=77a=quad('exp(-x/2).*sin(x+pi/6)',0,3*pi)a=0.9008

5/30/202427把matlab工作空間中一些有用的數(shù)據(jù)長久保存下來的方法是生成mat數(shù)據(jù)文件。

save——將工作空間中所有的變量存到matlab.mat文件中。數(shù)據(jù)的保存與獲取默認(rèn)文件名5/30/202428〔1〕savefilenamevariables將變量列表variables所列出的變量保存到磁盤文件filename中Variables所表示的變量列表中，不能用逗號(hào)，各個(gè)不同的變量之間只能用空格來分隔。未列出variables時(shí)，表示將當(dāng)前工作空間中所有變量都保持到磁盤文件中。(保存現(xiàn)場信息備用)缺省的磁盤文件擴(kuò)展名為“.mat”，數(shù)據(jù)存儲(chǔ)格式為二進(jìn)制?？梢允褂谩?”定義不同的存儲(chǔ)格式。5/30/202429savedata——將工作空間中所有的變量存到data.mat文件中。savedataab——將工作空間中a和b變量存到data.mat文件中。

下次運(yùn)行matlab時(shí)即可用load指令調(diào)用已生成的mat文件。5/30/202430〔2〕loadfilenamevariables將以前用save命令保存的變量