第2章 matlab

第二章 MATLAB绘图与插值拟合

§2.2 二维作图

1. 数值图

例2.2.3: t=0:0.1:2*pi; plotyy(t, sin(t),t,0.01*cos(t)) 例2.2.4

%li2_2_4.m n=1:100;

y1=(1+1./n).^n; y=(1+1./n).^(n+1); plot(n,y1,'.',n,y,'.r') 例2.2.5

function [x1,y1]= li2_2_5(n) %以li2_2_5.m文件名存盘 t1=1;x(1)=1;t2=2;y(1)=2; for k=1:n

x(k+1)=sqrt(x(k)*y(k)); y(k+1)=(x(k)+y(k))/2; end

x1=x(k+1);

y1=y(k+1); t=1:n+1;

plot(t,x,'-r',t,y,'.b') 例2.2.6

% li2_2_6.m

t=0:0.1:2*pi; y1= sin(t); y2=cos(t); y3=sin(t)*cos(t); plot(t, y1,'-.g', t,y2,':r',t,y3,'xb') 例2.2.7

%li2_2_7.m

x=-pi:pi/10:pi; y=tan(sin(x))-sin(tan(x));

52

plot(x, y,'--ro','LineWidth',2,'MarkerEdgeColor','k',… 'MarkerFaceColor','g',' MarkerSize',10)

例2.2.8 theta=0:0.1:8*pi; polar(theta,cos(5*theta)+1/3) 例2.2.9

%li2_2_9.m

x=-2:0.2:2;y=sin(x); subplot(2,2,1),plot(x,y) subplot(2,2,3),hist(y,5) subplot(2,2,4),stem(x,y);

subplot(2,2,2),rose(x,y); %stem( )和rose也是绘制二维曲线的函数。 Subplot(4,4,11), fill(x,y,'r'); Subplot(4,4,12), feather(x,y); Subplot(4,4,15), plot(x,y) Subplot(4,4,16), stairs(x,y);

例2.2.10 %li2_2_10.m subplot(2,2,1),

ezplot('5*x^4+2*x^2+1')

subplot(2,2,2),ezplot('sin(x)*exp(-x^2)') subplot(2,2,3),

ezplot('cos(x)/cos(2*x)') subplot(2,2,4),

ezplot('1+36*x/(x+3)^2')

§2.3 三维图形的绘制

例2.3.1

53

%li2_3_1.m

t=0:pi/50:10*pi x=sin(t);y=cos(t);z=t;h=plot3(x,y,z,'g-')

set(h,'LineWidth',4*get(h,'LineWidth')) 例2.3.2 %li2_3_2.m

[x,y]=meshgrid(-3:0.1:3,-3:0.1:3);%生成平面网格 z=x.^2-y.^2;

mesh(x,y,z) 例2.3.3

[x,y]=meshgrid(-3:0.1:3,-2:0.1:2); %li2_3_3.m z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); axis([-3,3,-2,2,-0.7,1.5])%设定坐标系范围 mesh(x,y,z)

§2.4 车灯光源投影区域的绘制（CUMCM 2002 A案例）

% li2_4_1.m p=0.03;x=25.0216;

for y1=-0.002:0.0004:0.002

y0=(-0.036:0.001:0.036)'*ones(1,73); z0=ones(73,1)*(-0.036:0.001:0.036); x0=(y0.^2+z0.^2)/(2*p);

xn=(p^3+4*x0*2*p.*x0+p*(-4*y1*y0+3*2*p*x0))./(2*(p^2+2*p*x0)); yn=(2*p*x0.*y0+p^2*(-y1+y0)+y1*(y0.^2-z0.^2))./(p^2+2*p*x0);

54

zn=(p^2+2*p*x0+2*y1*y0).*z0./(p^2+2*p*x0); y=y0+(yn-y0).*(x-x0)./(xn-x0); z=z0+(zn-z0).*(x-x0)./(xn-x0); plot(y,z,'b.') xlabel('y') ylabel('z') hold on end

§2.5 动画的绘制

例2.5.1 %li2_5_1.m clear

a=[-8/3 0 0;0 -10 10;0 28 -1]; y=[35 -10 -7]'; h=0.01;

p=plot3(y(1),y(2),y(3),'.','EraseMode','none','MarkerSize',5); %将擦除模式设置为none axis([0 50 -25 25 -25 25]) hold on for i=1:4000 a(1,3)=y(2); a(3,1)=-y(2); ydot=a*y; y=y+h*ydot;

set(p,'XData',y(1),'YData',y(2),'ZData',y(3))%设定图形目标的性质

55

drawnow%填充未完成的图形事件 i=i+1; end 例2.5.2 %li2_5_2.m clear

axis([0,2*pi,-1,1]) x=0:0.01:2*pi; plot(x,0) hold off hold on for j=1:30

t1=(j-1)*2*pi/30; t2=j*2*pi/30; t=t1:0.01:t2; plot(t,0,'.r') plot(t,sin(t),'.') end 例2.5.3 %li2_5_3.m clear

axis equal%创建一个用于显示该电影动画的坐标轴 m=moviein(16);

set(gca,'NextPlot','replacechildren') for j=1:16

plot(fft(eye(j+16))) m(:,j)=getframe; end movie(m,30)

56

%li2_7_1.m重金属污染地形图

A=xlsread('F:\\dimao.xls'); 读取F盘上地貌数据文件dimao.xls x=A(:,1); y=A(:,2); z=A(:,3);

scatter(x,y,5,z) %散点图 figure

[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4'); %插值 figure,surf(X,Y,Z,'EdgeColor','none') %三维曲面 title(‘重金属污染地形图’) grif off

%li2_7_2.m重金属浓度随地形分布图 clear all

A=xlsread('F:\\dimao.xls'); B=xlsread('F:\\ zjs.xls'); x=A(:,1); y=A(:,2);

z=B(:,1); %As数据，第i个重金属数据对应B矩阵中第i列 scatter(x,y,5,z) %散点图 figure

[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x),200)',linspace(min(y),max(y),200),'v4'); %插值 figure,surf(X,Y,Z,'EdgeColor','none') %三维曲面

title('As浓度随地域分布图') %As随对应重金属所在列做改动 grif off

62

§2.8拟合

例2.8.1编写如下程序

x= [0.5,1.0,1.5,2.0,2.5,3.0] %给出数据点的x值 y= [1.75，2.45,3.81,4.80,7.00,8.60]； %给出数据点的y值

p=polyfit(x,y,2)； %求出2阶拟合多项式的系数 poly2str(p,＇x＇) %求出2次拟和多项式表达式

＞＞ x=-2*pi:0.2:2*pi；y=1/5*cosh(x./3)； %将自变量和函数离散化 ＞＞ p=polyfit(x,y,2) %求出拟合多项式的系数 P=

0.0151 -0.0003 -.1836

＞＞ pk=poly2str(p,＇x＇) %求出拟合多项式表达式 pk=

0.015071x^2-0.00029557x+0.18365

＞＞ ezpolt 1/5*cosh(x/3)，hold %作悬链线函数曲线并保留 ＞＞ ezpot o.o15071*x^2-0.00029557*x+0.18365 %作拟合多项式函数曲线 ＞＞grid,legend(＇悬链线函数线＇，＇拟合悬链线的二次多项式曲线＇) 例2.8.3

function [a,b]=li2_8_3(c1,c2); %存盘文件名li2_8_3.m

%做非线性曲线拟合，参数传递时取c1=2;c2=1; x=1:5;

y=c1*exp(c2*x)+20*normrnd(0,1,1,5);%在指数曲线上加上随机误差，产生5个散点 [a,b]=lsqcurvefit(@nh,[c1,c2],x,y);% 非线性曲线拟合指令调用nh.m函数 x1=1:0.01:5; y2=c1*exp(c2*x1); y3=a(1)*exp(a(2)*x1);

63

plot(x1,y2,'b-',x,y,'ro',x1,y3,'g-'); %作散点图、原指数（y?c1eEnd

c2x）图、拟合图

function f=nh(a,x) %x数据通过主函数传递，a是由主函数传递过来的c1,c2迭代初值 %子函数程序，存盘文件名nh.m f=a(1)*exp(a(2).*x); end

例2.8.4 li2_8_4a.m.

function f=fun1(a,x) %li2_8_4a.m.

f=exp(a(1)*x+a(2)); % 用最小乘二拟合求上述函数中待定常数， %以及检验拟合效果的图形绘制程序 t=1790:10:1990;

x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76… 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4]; plot(t,x,’*’,t,x); a0=[0.001,1];

a=lsqcurvefit(‘fun1’,a0,t,x) ti=1790:10:2020; xi=fun1(a,ti); hold on plot(ti,xi); t1=2010; x1=fun1(a,t1) hold off

阻滞增长模型(logistic模型).

在MATLAB命令窗口键入：

64

Dsolve（’Dx=r*x*(1-x/c)’,’x(1790)=3.9’） ans=

c/(1+1/39*exp(-r*t)*exp(r)^1790*(10*c-39))

定义函数（2-8-3）的函数M文件. function f = fun3(a,t)

f=a(1)./(1+(a(1)/3.9-1)*exp(-(t-1790)*a(2)));

%li2_8_4b.m. x=1790:10:1990;

y=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76... 92 106.5 123.2 131.7 150.7 179.3 204 206.5 251.4]; plot(x,y,’*’,x,y); a0=[0.001,1];

a=lsqcurvefit(‘fun3’,a0,x,y) xi=1790:5:2020;yi=fun3(a,xi); hold on plot(xi,yi); x1=2010; y1=fun3(a,x1)

先编写一个函数M文件，以定义优化问题（2-8-4)式中的目标函数。 function f=fun5(a) n=16;w=2;

x=1790:10:1990;x1=x(1:n);x2=x(n+1:21);

y=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76... 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4]; y1=y(1:n);y2=y(n+1:21);

f=[fun3(a,x1),w*fun3(a,x2)]-[y1,w*y2];

主程序： %li2_8_4c.m.

65

t=1790:10:1990;

x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76... 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4]; plot(t,x,'*',t,x); a0=[300,0.03]; a=lsqnonlin('fun5',a0) ti=1790:5:1990; xi=fun3(a,ti); hold on; plot(ti,xi,t,x,'*'); x1=fun3(a,2010) hold off

例2.8.5

%Li2_8_5.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是向量. sd=c.*sd; s=[s,sum(sd)]; end

plot(z1,s,'r-'),hold on,xlabel('月份'),ylabel('平均气温')

66

§2.9黄河小浪底调水调沙问题

clc,clear all % Li2_9_1a.m

Load data3.txt%把表中5.8中的日期和时间数据行删除，余下的数据保存在纯文本文件 Liu=data3([1:3,:]);liu=liu’;liu=liu(:);%提出水流量并按照顺序变成列向量 sha=data3([2:4,:]);sha=sha’;sha=sha(:);%提出含沙量并按照顺序变成列向量 y=sha.*liu;y=y’;%计算排沙量，并变成行向量 i=1:24； t=(12*i-4)*3600; t1=t(1);t2=t(end);

pp=csape(t,p);%进行三次样条差值

xsh=pp.coefs%求得差值多项式的系数矩阵，每一行是一个区间上多项式的系数 TL=quadl（@（tt）ppval（pp，tt），t1，t2）

% Li2_9_1b.m clc,clear

load data3.txt%把表5.8中的日期和时间数据行删除，余下的数据保存在纯文本文件中 liu=data3([1,3],:);liu=liu';liu=liu(:);%提出水流量并按照顺序变成列向量 sha=data3([2,4],:);sha=sha';sha=sha(:);%提出含沙量并按照顺序变成列向量 y=sha.*liu;y=y';%计算排沙量，这是行向量 subplot(1,2,1),plot(liu(1:11),y(1:11),'*')

计算的Matlab程序如下： % Li2_9_1c.m clc,clear

load data3.txt%把表5.8中的日期和时间数据行删除，余下的数据保存在纯文本文件中 liu=data3([1:3,:]);liu=liu’;liu=liu(:);%提出水流量并按照顺序变成列向量

67

sha=data3([2:4,:]);sha=sha’;sha=sha(:);%提出含沙量并按照顺序变成列向量 y=sha.*liu;y=y’;%计算排沙量，这是列向量 format long e

%以下是第一阶段的拟合 for j=1:2

nihei{j}=polyfit(liu(1:11),y(1:11),j);%拟合多项式，系数排列从高次幂到低次幂 yhat1{j}=polyval(nihe1{j},liu(1:11));%求预测值 %以下求误差平方和与剩余标准差

cha1(j)=sum((y(1:11)-yhat1{j}).^2);rmsel(j)=sqrt(cha1(j)/(10-j)); end

nihe1{:}%显示细胞数组的所有元素 rmsel

%以下是第二阶段的拟合 for j=1:3

nihe2{j}=polyfit(liu(12:24),y(12:24),(j+1));%这里使用细胞数组 yhat{2}=polyval(nihe2{j},liu(12:24));

rmse2(j)=sqrt(sum((y(12:24)-yhat{j}).^2)/(11-j));%求剩余标准差 end

nihe{:}%显示细胞数组的所有元素 rmse2

format%恢复默认的短小数的现实形式

68

本文来源：https://www.bwwdw.com/article/x5vv.html