数学实验与 Matlab 321插值方法与多项式拟合的概念 322用 Matlab作插值和拟合 323用鼠标选节点观察插值、拟合的效果 324观察程序说明 ZXV 【 clf,a=-l;b=1;n=100 %用内联函数 inline命令定义函数, %在后面可直接用于函数g的计算,要改变函数做实验,可按此格式重新定义g g=inline('x2-x4): xX=linspace(a, b, n); gx(i)=gxx(i);%前面已经用 inline命令定义了g,可以这样用g计算函 数值 ymin=min(gx)*0.8;ymax=max(gx)*1.2;%分四个界面画图g的图形,以便于结 果比较 subplot(2, 2, 1) plot(xx,gx,-2),grid, hold on,axis([ a b ymin ymax]), title(近邻插值') ubplot(2, 2, 2) plot(xx,gx,-), grid. hold on,axis( a b ymin ymax]), title线性插值) subplot(2, 2, 3) plot(xx,gx,-2),grid, hold on,axis([ a b ymin ymax]), title('样条插值") subplot(2, 2, 4), plot(xx, gx, ' -- d, hold on,axis( a b ymin ymax],ttle多项式拟合) 用鼠标在屏幕上选点[x,y, button= input(n),可套用下面程序的格
数学实验与 Matlab 16 16 3.2.1 插值方法与多项式拟合的概念 3.2.2 用 Matlab 作插值和拟合 3.2.3 用鼠标选节点 观察插值、拟合的效果 3.2.4 观察程序说明 zxy3_1.m 【 clf,a=-1;b=1;n=100; %用内联函数 inline 命令定义函数， %在后面可直接用于函数 g 的计算 ,要改变函数做实验，可按此格式重新定义 g g=inline('x^2-x^4'); xx=linspace(a,b,n); for i=1:n gx(i)=gxx(i)); % 前面已经用 inline 命令定义了 g， 可以这样用 g 计算函 数 值 end ymin=min(gx)*0.8;ymax=max(gx)*1.2; %分四个界面画图 g 的图形，以便于结 果比较 subplot(2,2,1), plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('近邻插值') subplot(2,2,2), plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('线性插值') subplot(2,2,3), plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('样条插值') subplot(2,2,4),plot(xx,gx,'--'), grid,hold on,axis([a b ymin ymax]),title('多项式拟合') %用鼠标在屏幕上选点[x,y,button] = ginput(n)，可套用下面程序的格 式 button=1;
数学实验与 Matlab xl=a]; yl=[gx(1); while button==l [xi, yi, button]=ginput (1); subplot(2, 2, 1), h=plot(xi,yi,'ro') %在4个图形窗 口画点 subplot(2, 2, 2),h=plot(xi,yi, ro) subplot(2, 2, 3), h=plot(xi, yi, 'ro) subplot(2, 2, 4),h=plot(xi,yi, ro) I=[xi, x1; yl=yi,y1I %将选 的点存于向量x1,yl x1=[b, x1]; yl=[gx(n),y1]: XX=linspace(a, b, n) %定义自变量xx 计算不同的插值函数:x1,yl为节点,xx为输入自变量 nearest=interpI(xl, y 1, xx, 'nearest) linear=interpl(xl, yl, xx, 'linear) spline=interp(xl,y l, xX, ' spline) %多项式拟合指令[p,]= polyfit(x,y,n),n为拟合多项式次数,x,y为被拟合 数据 p为拟合多项式的系数,s是用来做误差估计和预测的数据结构 Ip, c]=polyfit(xl,y 1, 4) polyfit= polyol(p,xx);%用 polyol(p,x)计算系数为p的多项式在标量或 向量x处的值 subplot(2, 2, 1),h=plot(xx, nearest, 'r-); set(, ' linewidth, 2) %i subplot(2, 2, 2),h=plot(xx, linear, ' r-); set(h, 'linewidth, 2); subplot(2, 2, 3),h=plot(xx, spline, 'r-); set(h, 'linewidth, 2) subplot(2, 2, 4),h=plot(xx,ypoly fit, r-) set(h, ' linewidth, 2) 33应用、思考和练习
数学实验与 Matlab 17 17 x1=[a];y1=[gx(1)]; while button==1 [xi,yi,button]=ginput(1); subplot(2,2,1),h=plot(xi,yi,'ro') %在 4 个图形窗 口画点 subplot(2,2,2),h=plot(xi,yi,'ro') subplot(2,2,3),h=plot(xi,yi,'ro') subplot(2,2,4),h=plot(xi,yi,'ro') x1=[xi,x1];y1=[yi,y1]; %将 选 的点存于向量 x1,y1 end x1=[b,x1];y1=[gx(n),y1]; xx=linspace(a,b,n); %定义自变量 xx %计算不同的插值函数：x1,y1 为节点，xx 为输入自变量 ynearest=interp1(x1,y1,xx,'nearest'); ylinear=interp1(x1,y1,xx,'linear'); yspline=interp1(x1,y1,xx,'spline'); %多项式拟合指令[p,s] = polyfit(x,y,n)，n 为拟合多项式次数，x,y 为被拟合 数据， %p 为拟合多项式的系数，s 是用来做误差 估计和预测的数据结构。 [p,c]=polyfit(x1,y1,4); ypolyfit=polyval(p,xx); %用 polyval(p,x)计算系数为 p 的多项式在标量或 向 量 x 处的值 subplot(2,2,1),h=plot(xx,ynearest,'r-');set(h,'linewidth',2) %画 图 subplot(2,2,2),h=plot(xx,ylinear,'r-');set(h,'linewidth',2); subplot(2,2,3),h=plot(xx,yspline,'r-');set(h,'linewidth',2) subplot(2,2,4),h=plot(xx,ypolyfit,'r-');set(h,'linewidth',2) 】 3.3 应用、思考和练习
数学实验与 Matlab 331若干函数的插值和拟合练习 332几个应用问题 1.机床加工和水深流速问题 2.内燃机转角与升程的关系 3.随高度变化的大气压强 4.交通事故的调查 334多元函数的插值 1.矩形温箱的温度 2.航行区域的警示线 335 Fourier级数和周期函数的经验公式 zxy 3 2m n=10;,m=3;x=1:1:12; y=[3.13.86.912.716.820.524.525.922.016.110.75.4] z=pi/6)*; plot(z(1: 12),y(1: 12),'o); hold on %计算数据矩阵 X(1,)=[1cos(k*z(1))sin(k*z(1))] c=inv(X*X)*X"*y(l:n)’,%求解。 21=linspace (0, 2*pi, 30): s=[] %开始计算F一级数的部分和。 fori=1:30; sd=[1co(k*z1()sin(k*Z1(1));%注意k是向量
数学实验与 Matlab 18 18 3.3.1 若干函数的插值和拟合练习 3.3.2 几个应用问题 1. 机床加工和水深流速问题 2. 内燃机转角与升程的关系 3. 随高度变化的大气压强 4. 交通事故的调查 3.3.4 多元函数的插值 1. 矩形温箱的温度 2. 航行区域的警示线 3.3.5 Fourier级数和周期函数的经验公式 zxy3_2.m 【 clf,clear, n=10;m=3;x=1:1:12; y=[3.1 3.8 6.9 12.7 16.8 20.5 24.5 25.9 22.0 16.1 10.7 5.4]; z=(pi/6)*x;plot(z(1:12),y(1:12),'o');hold on k=1:m; %计算数据矩阵。 for i=1:n X(i,:)=[1 cos(k*z(i)) sin(k*z(i))]; end c=inv(X'*X)*X'*y(1:n)', %求 解 。 z1=linspace(0,2*pi,30); s=[]; %开 始计 算 F－级数的部分和。 for i=1:30; sd=[1 cos(k*z1(i)) sin(k*z1(i))]'; %注 意k是向量
数学实验与 Matlab d; s-[s, sum(sd); plot(z1,s,'r-), hold on, xlabel(月份), ylabel(平均气温 336由实验到理论:从开普勒定律和牛顿万有引力定律 实验4.微分、积分和微分方程 41.知识要点和背景:微积分学基本定理 42实验与观察():数值微积分 42.1实验:积分定义、微分方程和微积分基本定理的联系 ◆步骤1 t n=4 n1=n+ 1; xIlinspace(0, pi, nI ) fmyfun2 2(x): [0 n; x; f]) ◆步骤 y(1)=0;y(2)=y(1)+f(1)*(x(2)-x(1)); P intial=[x(1),y(1), P final=[x(2),y(2) 步骤3 【y(3)=(2)+f2)+(x(3)x(2) P intial=[x(2),y(3), P final=[x(3),y(3)
数学实验与 Matlab 19 19 sd=c.*sd; s=[s,sum(sd)]; end plot(z1,s,'r-'),hold on, xlabel('月份'),ylabel('平均气温 ') 】 3.3.6 由实验到理论:从开普勒定律和牛顿万有引力定律 实验４. 微分、积分和微分方程 4.1. 知识要点和背景：微积分学基本定理 4.2 实验与观察（Ⅰ）：数值微积分 4.2.1 实验：积分定义、微分方程和微积分基本定理的联系 ◆步骤1. 【 n=4;n1=n+1;x=linspace(0,pi,n1);f=myfun2_2(x);[0:n;x;f] 】 ◆步骤2. 【 y(1)=0; y(2)=y(1)+f(1) * (x(2)-x(1)); P_intial=[x(1),y(1)],P_final=[x(2),y(2)], 】 ◆ 步骤3. 【 y(3)=y(2)+f(2)* (x(3)-x(2)); P_intial=[x(2),y(3)], P_final=[x(3),y(3)] 】 ◆ 步骤4
数学实验与 Matlab Dy=y(3)-y(2)】 【Ds=f(2)*(x(3)-x(2))】 ◆步骤5 422求解数值积分的 Matlab积分命令 ◆观察 cumsum指令的功能。 【 ◆观察:用累积和命令 unum计算积分。 x=linspace(0, pi, 50); y=sin(x); T=cumsum(y)*pi/(50-1): 1=T(50) ◆观察梯形公式计算的效果。 x=linspace(0, pi, 50); y=sin(x): T=trapz(x, y) x ◆观察辛普森积分公式的计算效果 11=quad(sin, 0, pi), 12=quad(sin, 0, pi), 1 观察:用解常微分方程的方法计算分m( r y0=0: [x, y]=ode23('myfun2 2, [0, pi],y0); S-length(y); y(s) 1 【 y0=0; [x, y]=ode45('myfun2_ 2, [0, pily); s=length(y); y(s) I 423观察程序及其说明 观察方程的解曲线族,图4.1) clf, clear, a=0, b-pi, c-0: d=2. 2 n=20; [X,Y=meshgrid (linspace(a, b, n), linspace(c, d, n));
数学实验与 Matlab 20 20 【 Dy=y(3)-y(2) 】 【 DS=f(2)* (x(3)-x(2)) 】 ◆步骤5 . 4.2.2 求解数值积分的 Matlab 积分命令 ◆观察cumsum指令的功能。 【 x=[1 1 1 1 1]；I=cumsum(x) 】 ◆观察：用累积和命令cumsum计算积分。 【 x=linspace(0,pi,50); y=sin(x); T=cumsum(y)* pi/(50-1);I=T(50) 】 ◆观察梯形公式计算的效果。 【 x=linspace(0,pi,50); y=sin(x);T=trapz(x,y) 】 ◆ 观察辛普森积分公式的计算效果。 【 I1=quad('sin',0,pi), I2=quad8('sin',0,pi), 】 ◆ 观察：用解常微分方程的方法计算积分 0 sin(x)dx 。 【 y0=0;[x,y]=ode23('myfun2_2',[0,pi],y0);s=length(y);y(s) 】 【 y0=0;[x,y]=ode45('myfun2_2',[0,pi],y0);s=length(y);y(s) 】 4.2.3 观察程序及其说明 zxy4_1.m （观察方程的解曲线族，图4.1） 【 clf,clear,a=0;b=pi;c=0;d=2.2;n=20; [X,Y]=meshgrid(linspace(a,b,n),linspace(c,d,n));