数学实验与Matlab
式
16
button=1;
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应用、思考和练习
16
数学实验与Matlab
17
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是向量。
17
数学实验与Matlab
sd=c.*sd; s=[s,sum(sd)]; end
plot(z1,s,'r-'),hold on, xlabel('月份'),ylabel('平均气温') 】
18
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。
【 Dy=y(3)-y(2) 】
【 DS=f(2)*(x(3)-x(2)) 】
18
◆步骤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), 】
◆ 观察：用解常微分方程的方法计算积分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) 】
?0?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)); z=X.*Y; Fx=cos(atan(sin(X)));Fy=sin(atan(sin(X))); quiver(X,Y,Fx,Fy,0.5),hold on,axis([a,b,c,d]) [x,y]=ode45('zxy4_1f',[0,pi],0); [x1,y1]=ode45('zxy4_1f',[0,pi],0.2); [x2,y2]=ode45('zxy4_1f',[0,pi],0.4); [x3,y3]=ode45('zxy4_1f',[0,pi],0.6);
plot(x,y,'r.-',x1,y1,'r.-',x2,y2,'r.-',x3,y3,'r.-'), xlabel('x') ,ylabel('y') 】 zxy4_1f.m
19
【 function dy=zxy4_1f(x,y)
dy=sin(x); 】 . 观察 程序zxy4_2.m （图4.1～4.3）
【 clf, a=0;b=4*pi;n=31;
x=linspace(a,b,n); y=zeros(1,n); y(1)=0; s(1)=0;
for i=1:n-1 dy=myfun2_2(x(i)); y(i+1)=y(i)+dy*(x(i+1)-x(i)); hh(i)=myfun2_2(x(i));
s(i+1)=s(i)+hh(i)*(x(i+1)-x(i)); end
a1=0.9*a;b1=1.1*b; % 设置坐标范围。 ymin=min(y');ymax=max(y'); ymin1=ymin*0.9;ymax1=ymax*1.1;
subplot(2,1,1), %在第一幅图中画垂线,和原函数。
fplot('1-cos(x)',[0,pi],'r:');axis([a1 b1 ymin1 ymax1]),hold on, plot([x(1) x(1)],[ymin ymax]);
subplot(2,1,2), %在第二幅图中画被积函数图象。
fplot('myfun2_2',[a b],'-k');hold on,axis([a1 b1 ymin1 ymax1])
for i=1:n-1 %在第I个子区间上计算。
subplot(2,1,1),
plot([x(i+1) x(i+1)],[ymin ymax]); %在各分点画竖线。 plot([x(i) x(i+1)],[y(i),y(i)]); %画水平线。 h=plot([x(i+1) x(i+1)],[y(i),y(i+1)]); %画表示增量的垂线。 set(h,'linewidth',3,'color','b'); %设置线宽和颜色。
h1=plot([x(i) x(i+1)],[y(i) y(i+1)],'.-'); %画折线，设置图形属性 。 set(h1,'linewidth',2,'markersize',18); subplot(2,1,2),
xx=[x(i) x(i) x(i+1) x(i+1)]; yy=[0 hh(i) hh(i) 0]; %计算矩形顶
点坐标。
patch(xx,yy,[0.7 0.7 0.7]); %在第二幅图中画矩形块并填充颜色。
plot([x(i+1) x(i+1)],[ymin ymax]);
plot([x(1) x(1)],[s(i),s(i+1)],'r','linewidth',5); %在y轴上画面积增
20