MATLAB编程基础第6讲-插值、拟合与初值常微分方程的求解

MATLAB编程根底

之插值、拟合与初值常微分方程的求解第六讲梁丙臣5/30/20241多项式的拟合需求：实验数据总结为规律曲线等等最小二乘法：所用曲线限定为多项式，在数据点上拟合值和函数值有最小误差平方和。P=polyfit(x,y,n)：把自变量x和函数值y拟合成n阶多项式，向量为p。5/30/20242例3-35x=(0:0.1:1)';y=[-0.41.93.26.27.17.37.79.69.59.312]';%给出一组11个点数据y2=polyfit(x,y,2)%计算2阶拟合的多项式向量

x1=0:0.01:1;f2=polyval(y2,x1);%2阶拟合曲线在各点的函数值y10=polyfit(x,y,10)%计算10阶拟合的多项式向量

f10=polyval(y10,x1);%10阶拟合曲线在各点的函数值plot(x,y,'o',x1,f2,':',x1,f10,'k')5/30/20243例3-36x=(0:0.1:2.5)';%给出一组数据，为误差函数的一个区间y=erf(x);p=polyfit(x,y,6)%计算该区间内6阶拟合多项式的向量

x1=(0:0.1:5)';%将区间增长一倍y1=erf(x1);%计算误差函数在新区间内的函数值f=polyval(p,x1);%计算6阶拟合曲线在新区间内的取值plot(x1,y1,'o',x1,f,'-')%绘图，比较前后区间内的曲线拟合效果，如图3-7axis([0502])5/30/20244例：拟合以下数据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例：根据经验公式y=a+bx2，拟合如下数据：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多项式的插值插值的定义——是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。Matlab提供了一维、二维、三次样条等许多插值选择5/30/202471.一维插值yi=interp1(x,Y,xi,method)：求同维数据x和Y，运用method指定的方法计算插值点xi处的数值yi。method有如下四种：nearest最近点差值，取与最近的数据点的值linear线性差值，直线连接数据点，插值点值位于直线上spline样条插值，用三次样条曲线通过数据点，根据曲线进行差值cubic立方差值，用三次曲线拟合并通过数据点5/30/20248例3-37x=0:10;%给出基准数据y=sin(x);xi=0:.25:10;%在两个基准数据点间插入3个点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--')%绘图，并标注各曲线代表的插值方法，如图3-8legend('原始数据','最近点插值','线性插值','样条插值','立方插值')5/30/202492.二维插值ZI=interp2(X,Y,Z,XI,YI,method)：同维数据X,Y和Z，运用method指定的方法计算插值点(XI,YI)处的数值ZI。5/30/202410例3-38axis([-33-33-520])%确定坐标轴的范围[X,Y]=meshgrid(-3:.5:3);%生成网格Z=peaks(X,Y);%计算peaks函数值figure(1)mesh(X,Y,Z)%绘制原始数据曲面图，如图3-9

5/30/202411[XI,YI]=meshgrid(-3:.125:3);%确定插值点ZI1=interp2(X,Y,Z,XI,YI,'*nearest');%计算用nearest方法所得二维插值figure(2)mesh(XI,YI,ZI1)%绘制最近点插值法曲面图，如图3-10ZI2=interp2(X,Y,Z,XI,YI,'*linear');%计算用linear法所得二维插值figure(3)mesh(XI,YI,ZI2)%绘制双线性插值法曲面图，如图3-115/30/202412ZI3=interp2(X,Y,Z,XI,YI,'*spline');%计算用spline法所得二维插值figure(4)mesh(XI,YI,ZI3)%绘制双样条插值法曲面图，如图3-12ZI4=interp2(X,Y,Z,XI,YI,'*cubic');%计算用cubic法所得二维插值figure(5)mesh(XI,YI,ZI4)%绘制双立方插值法曲面图，如图3-135/30/2024133.9常微分方程初值问题的数值解法龙格－库塔法简介龙格－库塔法的实现基于龙格－库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：[t,y]=ode23('fname',tspan,y0)[t,y]=ode45('fname',tspan,y0)其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。5/30/202414一阶常微分方程的初值问题求解例3-39，求微分方程组在区间[012]内的数值解，且满足初始条件5/30/202415例3-39%首先编写方程函数，名为rigid.mfunctiondy=rigid(t,y)dy=zeros(3,1);%生成3×1的矩阵dy(1)=y(2)*y(3);%第一个元素是dy(2)=-y(1)*y(3);%第二个元素是dy(3)=-0.51*y(1)*y(2);%第三个元素是

5/30/202416例3-39〔续前〕tspan=[012];y0=[011];%在同一目录下，计算方程数值解。输入时间区间和初始条件[t,Y]=ode45('rigid',tspan,y0);%采用ode45算法求解方程，options为默认值plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')%绘制计算结果并标注，如图3-14legend('Y1','Y2','Y3')5/30/202417高阶微分方程的初值问题求解高阶微分方程：把高阶微分方程转换为一阶方程组，初始条件也做相应的替换。通常令：上式改写为：5/30/202418例3-40，求方程初始条件在时间区间[012]内的数值解把高阶微分转化低阶微分形式，令：原方程变为：初始条件：5/30/202419例3-40训练任务：请编制名为odet3.m方程函数tspan=[012];%在同一目录下，计算方程数值解。给出时间区间和初始值y0=[-110];[t,Y]=ode45('odet3',tspan,y0);%采用ode45算法plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')%绘制计算结果并标注，如图3-15legend('Y1','Y2','Y3')5/30/202420%首先编写方程函数，名为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例设有初值问题，试求其数值解，并与精确解相比较(精确解为y(t)=)。(1)建立函数文件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);%求数值解y1=sqrt(t+1)+1;%求精确解t'y'y1'y为数值解，y1为精确值，显然两者近似。5/30/202422例：x+(x2-1)x+x=0为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式

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)给定区间、初始值;求解微分方程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训练1a=[34-7-12;5-742;108-5;-65-210];b=[4-39-8]';c=a\b5/30/202426训练2训练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工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。

save——将工作空间中所有的变量存到matlab.mat文件中。数据的保存与获取默认文件名5/30/202428〔1〕savefilenamevariables将变量列表variables所列出的变量保存到磁盘文件filename中Variables所表示的变量列表中，不能用逗号，各个不同的变量之间只能用空格来分隔。未列出variables时，表示将当前工作空间中所有变量都保持到磁盘文件中。(保存现场信息备用)缺省的磁盘文件扩展名为“.mat”，数据存储格式为二进制。可以使用“-”定义不同的存储格式。5/30/202429savedata——将工作空间中所有的变量存到data.mat文件中。savedataab——将工作空间中a和b变量存到data.mat文件中。

下次运行matlab时即可用load指令调用已生成的mat文件。5/30/202430〔2〕loadfilenamevariables将以前用save命令保存的变量