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

2020/2/101MATLAB编程基础之插值、拟合与初值常微分方程的求解第六讲梁丙臣2020/2/102•3.8.3多项式的拟合需求：实验数据总结为规律曲线等等最小二乘法：所用曲线限定为多项式，在数据点上拟合值和函数值有最小误差平方和。P=polyfit(x,y,n)：把自变量x和函数值y拟合成n阶多项式，向量为p。2020/2/103•例3-35•x=(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')2020/2/104•例3-36•x=(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-7•axis([0502])2020/2/105例：拟合以下数据x0.51.01.52.02.53.0y1.752.453.814.808.008.60x=[0.51.01.52.02.53.0];y=[1.752.453.814.808.008.60];a=polyfit(x,y,2)a=0.49001.25010.8560x1=[0.5:0.05:3.0];y1=a(3)+a(2)*x1+a(1)*x1.^2;plot(x,y,'*')holdonplot(x1,y1,'-r')0.511.522.5302468102020/2/106例：根据经验公式y=a+bx2，拟合如下数据：xi1925313844yi19.032.349.073.398.8x=[1925313844];y=[19.032.349.073.398.8];x1=x.^2;x1=[ones(5,1),x1'];ab=x1\y'ab=0.59370.0506x0=[19:0.2:44];y0=ab(1)+ab(2)*x0.^2;plot(x,y,'o',x0,y0,'-r')152025303540450204060801002020/2/1073.8.4多项式的插值插值的定义——是对某些集合给定的数据点之间函数的估值方法。•当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。•Matlab提供了一维、二维、三次样条等许多插值选择2020/2/108•1.一维插值yi=interp1(x,Y,xi,method)：求已知同维数据x和Y，运用method指定的方法计算插值点xi处的数值yi。method有如下四种：nearest最近点差值，取与最近的已知数据点的值linear线性差值，直线连接数据点，插值点值位于直线上spline样条插值，用三次样条曲线通过数据点，根据曲线进行差值cubic立方差值，用三次曲线拟合并通过数据点2020/2/109•例3-37•x=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-8•legend('原始数据','最近点插值','线性插值','样条插值','立方插值')2020/2/10102.二维插值ZI=interp2(X,Y,Z,XI,YI,method)：已知同维数据X,Y和Z，运用method指定的方法计算插值点(XI,YI)处的数值ZI。2020/2/1011•例3-38•axis([-33-33-520])•%确定坐标轴的范围•[X,Y]=meshgrid(-3:.5:3);•%生成网格•Z=peaks(X,Y);•%计算peaks函数值•figure(1)•mesh(X,Y,Z)•%绘制原始数据曲面图，如图3-9••2020/2/1012•[XI,YI]=meshgrid(-3:.125:3);•%确定插值点•ZI1=interp2(X,Y,Z,XI,YI,'*nearest');•%计算用nearest方法所得二维插值•figure(2)•mesh(XI,YI,ZI1)•%绘制最近点插值法曲面图，如图3-10•ZI2=interp2(X,Y,Z,XI,YI,'*linear');•%计算用linear法所得二维插值•figure(3)•mesh(XI,YI,ZI2)•%绘制双线性插值法曲面图，如图3-112020/2/1013ZI3=interp2(X,Y,Z,XI,YI,'*spline');•%计算用spline法所得二维插值•figure(4)•mesh(XI,YI,ZI3)•%绘制双样条插值法曲面图，如图3-12•ZI4=interp2(X,Y,Z,XI,YI,'*cubic');•%计算用cubic法所得二维插值•figure(5)•mesh(XI,YI,ZI4)•%绘制双立方插值法曲面图，如图3-132020/2/10143.9常微分方程初值问题的数值解法龙格－库塔法简介龙格－库塔法的实现基于龙格－库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为：[t,y]=ode23('fname',tspan,y0)[t,y]=ode45('fname',tspan,y0)其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。2020/2/1015•一阶常微分方程的初值问题求解00,,yxybxayxfy21331232151.0yyyyyyyyy•例3-39，求微分方程组在区间[012]内的数值解，且满足初始条件101000321yyy2020/2/1016•例3-39•%首先编写方程函数，名为rigid.m•functiondy=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);•%第三个元素是•2020/2/1017•例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-14•legend('Y1','Y2','Y3')2020/2/1018•高阶微分方程的初值问题求解0,,,,,xyyyyFn101000,0,0nnyyyyyy121,,,nnyyyyyyYxfYxfYxfyyyYnn,,,21211000210000nnyyyyyyY高阶微分方程：把高阶微分方程转换为一阶方程组，初始条件也做相应的替换。通常令：上式改写为：2020/2/1019•例3-40，求方程初始条件042212yyyy00,10,10yyy在时间区间[012]内的数值解把高阶微分转化低阶微分形式，令：yyyyyy321,,原方程变为：4222321133221yyyyyyyy初始条件：001010321yyy2020/2/1020•例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-15•legend('Y1','Y2','Y3')2020/2/1021•%首先编写方程函数，名为odet3.m•functiondy=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);2020/2/1022例设有初值问题，试求其数值解，并与精确解相比较(精确解为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为精确值，显然两者近似。2020/2/1023例：x+(x2-1)x+x=0为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式x1=x2x2=x2(1-x12)-x11.建立m文件2.解微分方程······2020/2/1024建立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))2020/2/1025-2.5-2-1.5-1-0.500.511.522.5-3-2-10123051015202530-3-2-101232020/2/1026训练181025695832475412743wzyxwzxwzyxwzyxa=[34-7-12;5-742;108-5;-65-210];b=[4-39-8]';c=a\b2020/2/1027训练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=77320sin(/6)xexpidxa=quad('exp(-x/2).*sin(x+pi/6)',0,3*pi