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

1、2020/7/7,1,MATLAB编程基础之,插值、拟合与初值常微分方程的求解,第六讲,梁丙臣,2020/7/7,2,3.8.3多项式的拟合 需求：实验数据总结为规律曲线等等 最小二乘法：所用曲线限定为多项式，在数据点上拟合值和函数值有最小误差平方和。 P=polyfit(x,y,n)：把自变量x和函数值y拟合成n阶多项式，向量为p。,2020/7/7,3,例3-35 x = (0:0.1:1); y = -0.4 1.9 3.2 6.2 7.1 7.3 7.7 9.6 9.5 9.3 12; % 给出一组11个点数据 y2 = polyfit(x,y,2) % 计算2阶拟合的多项式向量 x1

2、= 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/7/7,4,例3-36 x = (0: 0.1: 2.5); % 给出一组数据，为误差函数的一个区间 y = erf(x); p = polyfit(x,y,6) % 计算该区间内6阶拟合多项式的向量 x1 = (0: 0.1: 5); % 将区间增长一倍 y1 = erf(x

3、1); % 计算误差函数在新区间内的函数值 f = polyval(p,x1); % 计算6阶拟合曲线在新区间内的取值 plot(x1,y1,o,x1,f,-) % 绘图，比较前后区间内的曲线拟合效果，如图3-7 axis(0 5 0 2),2020/7/7,5,例：拟合以下数据,x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 8.00 8.60; a=polyfit(x,y,2) a = 0.4900 1.2501 0.8560 x1=0.5:0.05:3.0; y1=a(3)+a(2)*x1+a(1)*x1.2; plot(x,y,*) h

4、old on plot(x1,y1,-r),2020/7/7,6,例：根据经验公式y=a+bx2，拟合如下数据：,x=19 25 31 38 44; y=19.0 32.3 49.0 73.3 98.8; x1=x.2; x1=ones(5,1),x1; ab=x1y ab = 0.5937 0.0506 x0=19:0.2:44; y0=ab(1)+ab(2)*x0.2; plot(x,y,o,x0,y0,-r),2020/7/7,7,3.8.4多项式的插值 插值的定义是对某些集合给定的数据点之间函数的估值方法。 当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。 Matlab

5、提供了一维、二维、 三次样条等许多插值选择,2020/7/7,8,1.一维插值 yi=interp1(x,Y,xi,method)：求已知同维数据x和Y，运用method指定的方法计算插值点xi处的数值yi。 method有如下四种： nearest 最近点差值，取与最近的已知数据点的值 linear 线性差值，直线连接数据点，插值点值位于 直线上 spline 样条插值，用三次样条曲线通过数据点，根据曲线进行差值 cubic 立方差值，用三次曲线拟合并通过数据点,2020/7/7,9,例3-37 x = 0:10; % 给出已知基准数据 y = sin(x); xi = 0:.25:10;

6、% 在两个基准数据点间插入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/7/7,10,2.二维插值 ZI=interp2

7、(X,Y,Z,XI,YI,method)：已知同维数据X,Y和Z，运用method指定的方法计算插值点(XI,YI)处的数值ZI。,2020/7/7,11,例3-38 axis(-3 3 -3 3 -5 20) % 确定坐标轴的范围 X,Y = meshgrid(-3:.5:3); % 生成网格 Z = peaks(X,Y); % 计算peaks函数值 figure(1) mesh(X,Y,Z) % 绘制原始数据曲面图，如图3-9,2020/7/7,12,XI,YI = meshgrid(-3:.125:3); % 确定插值点 ZI1 = interp2(X,Y,Z,XI,YI,*neares

8、t); % 计算用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-11,2020/7/7,13,ZI3 = interp2(X,Y,Z,XI,YI,*spline); % 计算用spline法所得二维插值 figure(4) mesh(XI,YI,ZI3) % 绘制双样条插值法曲面图，如图3-12 ZI4 = interp

9、2(X,Y,Z,XI,YI,*cubic); % 计算用cubic法所得二维插值 figure(5) mesh(XI,YI,ZI4) % 绘制双立方插值法曲面图，如图3-13,2020/7/7,14,3.9常微分方程初值问题的数值解法 龙格库塔法简介 龙格库塔法的实现 基于龙格库塔法，MATLAB提供了求常微分方程数值解的函数，一般调用格式为： t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为t0,tf,表示求解区间。y0是初始状态列向量。t和y分

10、别给出时间向量和相应的状态向量。,2020/7/7,15,一阶常微分方程的初值问题求解,例3-39，求微分方程组,在区间0 12内的数值解，且满足初始条件,2020/7/7,16,例3-39 % 首先编写方程函数，名为rigid.m function dy = rigid(t,y) dy = zeros(3,1); % 生成31的矩阵 dy(1) = y(2) * y(3); % 第一个元素是 dy(2) = -y(1) * y(3); % 第二个元素是 dy(3) = -0.51 * y(1) * y(2); % 第三个元素是,2020/7/7,17,例3-39（续前） tspan = 0

11、12; y0 = 0 1 1; % 在同一目录下，计算方程数值解。输入时间区间和初始条件 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/7/7,18,高阶微分方程的初值问题求解,高阶微分方程：,把高阶微分方程转换为一阶方程组，初始条件也做相应的替换。通常令：,上式改写为：,2020/7/7,19,例3-40，求方程,初始条件,在时间区间0 12内的数值解,把高阶微分转

12、化低阶微分形式，令：,原方程变为：,初始条件：,2020/7/7,20,例3-40 训练任务：请编制名为odet3.m方程函数 tspan = 0 12; % 在同一目录下，计算方程数值解。给出时间区间和初始值 y0 = -1 1 0; t,Y = ode45(odet3,tspan,y0); % 采用ode45算法 plot(t,Y(:,1),-,t,Y(:,2),-.,t,Y(:,3),.) % 绘制计算结果并标注，如图3-15 legend(Y1,Y2,Y3),2020/7/7,21,% 首先编写方程函数，名为odet3.m function dy = odet3(t,y) dy = z

13、eros(3,1); dy(1) = y(2); dy(2) = y(3); dy(3) = y(1).(1/2)-2*y(3).2*(y(2)-4);,2020/7/7,22,例设有初值问题，试求其数值解，并与精确解相比较(精确解为y(t)=)。 (1) 建立函数文件funt.m。 function yp=funt(t,y) yp=(y2-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为精确值，显然两者近似。,20

14、20/7/7,23,例：x+(x2-1)x+x=0 为方便令x1=x，x2=x分别对x1,x2求一 阶导数，整理后写成一阶微分方程组 形式 x1=x2 x2=x2(1-x12)-x1 建立m文件 解微分方程,2020/7/7,24,建立m文件 function xdot=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=0 0.25; t,x=ode23(wf, t0, tf, x0) figure(1), plot(t,x)； figure(2),p

15、lot(x(:,1),x(:,2),2020/7/7,25,2020/7/7,26,训练1,a=3 4 -7 -12;5 -7 4 2;1 0 8 -5;-6 5 -2 10; b=4 -3 9 -8; c=ab,2020/7/7,27,训练2 训练3,算出 x3+2x2+x+1=0 的根,function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); S,n=quad(fesin,0,3*pi) S = 0.9008 n = 77,a=quad(exp(-x/2).*sin(x+pi/6),0,3*pi),0.9008,2020/7/7,28,把matlab工

16、作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。 save 将工作空间中所有的变量存到matlab.mat文件中。,数据的保存与获取,默认文件名,2020/7/7,29,（1）save filename variables 将变量列表variables所列出的变量保存到磁盘文件filename中 Variables所表示的变量列表中，不能用逗号，各个不同的变量之间只能用空格来分隔。 未列出variables时，表示将当前工作空间中所有变量都保持到磁盘文件中。(保存现场信息备用) 缺省的磁盘文件扩展名为“.mat”，数据存储格式为二进制。可以使用“-”定义不同的存储格式。,2020/7/7,30,save data将工作空间中所有的变量存到data.mat文件中。 save data a b 将工作空间中a和b变量存到data.mat文件中。 下次运行matlab时即可用load指令调用已生成的mat文件。,2020/7/7,31,（2）load filename variables 将以前用save命令