MATLAB应用实例分析例分析Word下载.docx
《MATLAB应用实例分析例分析Word下载.docx》由会员分享,可在线阅读,更多相关《MATLAB应用实例分析例分析Word下载.docx(11页珍藏版)》请在冰豆网上搜索。
i=0.5:
3.5;
p=1./(n*log(1+0.01*i));
T=log(r'
)*p;
plot(r,T)
xlabel('
r'
)%给x轴加标题
ylabel('
T'
)%给y轴加标题
q=ones(1,length(i));
text(7*q-0.2,[T(14,1:
5)+0.5,T(14,6)-0.1,T(14,7)-0.9],num2str(i'
))
40
350.530
25
20T115
1.51022.5533.50123456789r
图1
1
从图1中既可以看到T随r的变化规律,而且还能看到i的不同取值对T—r曲线的影响(图中的六条曲线分别代表i的不同取值)。
二、已知多项式求根
65432已知多项式为,求其根。
h,x,10x,31x,10x,116x,200x,96
对多项式求根问题,我们常用roots()函数。
>
h=roots([1-1031-10-116200-96])%中括号内为多项式系数由高阶到常数。
计算结果显示为(其中i为虚数单位):
h=
-2.0000
4.0000
3.0000
2.0000+0.0000i
2.0000-0.0000i
1.0000
如果已知多项式的根,求多项式,用poly()函数。
对上面得到的h的值求多项式,其MATLAB的表达形式及结果如下:
h=[-2.00004.00003.00002.0000+0.0000i2.0000-0.0000i1.0000];
c=poly(h)
计算结果显示为:
c=
1-1031-10-116200-96
三、方程组的求解
8x,x,6x,7.5,123,求解下面的方程组:
3x,5x,7x,4,123,4x,9x,2x,12123,
对于线性方程组求解,常用线性代数的方法,把方程组转化为矩阵进行计算。
1,x,abax,b,x,a\b
a=[816;
357;
492];
%建立系数矩阵
b=[7.5;
4;
12];
%建立常数项矩阵
x=a\b%求方程组的解
x=
1.2931
0.8972
-0.6236
四、数据拟合与二维绘图
在数学建模竞赛中,我们常会遇到这种数据表格问题,如果我们仅凭眼睛观察,很难看到其中的规律,也就更难写出有效的数学表达式从而建立数学模型。
因此可以利用MATLAB的拟合函数,即polyfit()
函数,并结合MATLAB的绘图功能(利用plot()函数),得到直观的表示。
2
例:
在化学反应中,为研究某化合物的浓度随时间的变化规律,测得一组数据如下表:
12345678T(分)
y46.48.08.49.289.59.79.86
910111213141516T(分)
y1010.210.3210.4210.510.5510.5810.6
MATLAB的表达形式如下:
t=[1:
16];
%数据输入
y=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];
plot(t,y,'
o'
)%画散点图
p=polyfit(t,y,2)%二次多项式拟合
holdon
xi=linspace(0,16,160);
%在[0,16]等间距取160个点
yi=polyval(p,xi);
%由拟合得到的多项式及xi,确定yi
plot(xi,yi)%画拟合曲线图
执行程序得到图2;
11
10
9
8
7
6
5
40246810121416
图2
显示的结果为
p=
0.04451.07114.3252-
2p的值表示二阶拟合得到的多项式为:
y=-0.0445t+1.0711t+4.3252
下面是用lsqcurvefit()函数,即最小二乘拟合方法的Matlab表达:
x0=[0.1,0.1,0.1];
zuixiao=inline('
x
(1)*t.^2+x
(2)*t+x(3)'
'
x'
t'
);
x=lsqcurvefit(zuixiao,x0,t,y)%利用最小二乘拟合
其显示的结果为:
-0.04451.07114.3252
可以看出其得到的结果与polyfit函数的结果相同。
这说明在多项式拟合问题上这两个函数的效果是相同的。
3
下面的一个例子将体现lsqcurvefit()函数的优势。
例2:
在物理学中,为研究某种材料应力与应变的关系,测得一组数据如下表:
925112516252125262531253625应力σ
0.110.160.350.480.610.710.85应变ε
如果假定应力与应变有如下关系(σ为应力值,ε为应变值):
ε=a+blnσ试计算a、b的值。
x=[925,1125,1625,2125,2625,3125,3625];
y=[0.11,0.16,0.35,0.48,0.61,0.71,0.85];
plot(x,y,'
)
[p,resid1]=polyfit(x,y,2)
xi=linspace(700,3700,3000);
yi=polyval(p,xi);
plot(xi,yi)
x0=[0.1,0.1];
fff=inline('
a
(1)+a
(2)*log(x)'
a'
[a,resid2]=lsqcurvefit(fff,x0,x,y)plot(xi,fff(a,xi),'
执行程序得到图3,图中蓝色曲线为利用polyfit()函数得到的曲线,红色曲线为利用lsqcurvefit()函
数得到的曲线;
0.90.80.70.60.50.40.30.20.10-0.15001000150020002500300035004000
p=
-0.00000.0004-0.2266resid1=
R:
[3x3double]
df:
4
normr:
0.0331
a=
-3.58100.5344
resid2=
0.0064
其中a的值代表利用lsqcurvefit()函数得到的关系为:
ε=-3.5810+0.5344σ
4
resid1、resid2分别代表运用polyfit()函数、lsqcurvefit()函数得到的残差。
可以看出利用lsqcurvefit()函数残差更小,即得到了更好的拟合效果。
在数学建模的实际问题中,如果问题的机理不明,我们只能采用polyfit()函数,即多项式拟合的
方法,以获得近似的数据描述函数;
但如果通过分析,可以得到一些机理,那么采用最小二乘的方法将
得到更好的效果,而且得到的拟合函数也更有意义。
五、隐函数的图形绘制
1plot()只能绘制显函数图形,对于形如的复杂隐函数,很难转化为,lny,ln(,1,y),x,sin(x),0y
显函数并利用plot()函数绘制图形,这时就可以用ezplot()函数直接绘制其曲线。
MATLAB的表达形式如下:
ezplot('
1/y-log(y)+log(-1+y)+x-sin(x)'
执行程序得到图5
x=sin(3t)cos(t),y=sin(3t)sin(t)1/y-log(y)+log(-1+y)+x-sin(x)=060.4
40.2
20
-0.2yy0
-0.4-2-0.6-4-0.8
-6-0.8-0.6-0.4-0.200.20.40.60.8-6-4-20246xx
图5图6如果是形如下面的参数方程,同样可以利用ezplot()函数绘制其曲x,sin3tcost,y,sin3tsint,t,(0,,)线。
MATLAB的表达形式如下:
sin(3*t)*cos(t)'
sin(3*t)*sin(t)'
[0,pi])
执行程序得到图6。
六、三维图形绘制
假设有一个时间向量t,对该向量进行下列运算则可以构成三个坐标值向量
x,sint,y,cost,z,t
对于上面的方程可以利用ezplot3()函数或plot3()函数绘制三维曲线。
这里仅列举ezplot3()函数的使用。
ezplot3('
sin(t)'
cos(t)'
[0,6*pi])
执行程序得到图7:
3绘制下述曲面:
z(r,,),rcos(3,),其中0,r,1,0,,,2,
nr=12;
nth=50;
r=linspace(0,1,nr);
theta=linspace(0,2*pi,nth);
[R,T]=meshgrid(r,theta)
x=cos(theta'
)*r;
5
y=sin(theta'
surf(x,y,R.^3.*cos(3*T))
执行程序得到图8。
x=sin(t),y=cos(t),z=t
15
10z
010.510.500-0.5-0.5-1-1yx
图7图8
除了surf()函数还有surfc()、surfl()、mesh()、waterfall()函数也用于曲面的绘制,具体效果如图9所示,可
以针对自己的需要选取适合的曲面绘制函数。
图9
MATLAB作图之三维绘图示例山体绘制
%三维绘图示例山体绘制
%mesh函数演示
x=1.0:
0.1:
2.0;
y=2.0:
3.0;
[X,Y]=meshgrid(x,y);
z=[5.115.135.145.135.095.044.984.934.894.854.85
5.395.495.515.465.325.144.944.744.594.494.48
5.615.775.815.715.515.234.904.594.364.214.19
5.735.925.975.865.625.274.884.514.234.054.03
5.745.925.975.867.625.274.884.514.214.044.02
5.635.795.846.7410.539.238.914.594.334.184.16
5.425.535.565.497.355.164.934.734.554.454.44
5.145.185.195.1711.125.054.974.904.844.814.80
4.484.804.794.824.874.945.025.105.165.195.20
4.564.454.434.494.644.845.065.285.455.555.56
4.364.194.164.254.474.765.095.415.665.815.83];
figure
(1)
%画网格图
mesh(X,Y,z);
6
colormap([010]);
x轴'
y轴'
zlabel('
z轴'
%画表面图
figure
(2)
surf(X,Y,z);
colormap([100]);
1212
1010
88z轴z轴66
4433221.81.82.52.51.61.61.41.41.21.22211y轴y轴x轴x轴
七、二项分布的使用
飞机成功起飞的概率问题:
由16架飞机组成的空军飞行中队要求做好立即起飞的准备,其中一架飞机不能立即起飞的概率为
20%,重新起飞需几分钟的时间,因此一架飞机立刻起飞的概率为0.80。
12架飞机能够成功起飞的概率
为多少,
这是一个概率中的二项分布问题,常用binopdf()函数。
h=binopdf(12,16,0.80)%二项分布函数的概率值
计算结果为
h=0.2001
另一方面,至少有14架飞机立刻成功起飞的概率为:
h=1-binocdf(13,16,0.80)%或h=sum(binopdf(14:
16,16,0.80)),其中binocdf()为二项分布的累积概率值,计算结果为h=0.3518。
在实际的数学建模竞赛中,仅罗列一个一个的数据是枯燥而又不直观的,很难吸引人们的注意,也不容
易打动评委们的心;
因此,结合数值计算结果,并合理的利用MATLAB的绘图功能会起到事半功倍的效
果。
下面的程序为运行结果的绘图(图10)表示:
n=1:
16;
h=binopdf(n,16,0.80);
plot([n;
n],[zeros(1,16);
h],'
k'
)%二维绘图函数
text(8-.7:
16-.7,h(8:
16)+.005,num2str(h(8:
16)'
3))%在图中进行注释函数axis([01700.27])%坐标轴取值范围函数
Numberofaircraftlaunchedontime'
)%给x轴加标题ylabel('
probability'
)%给y轴加标题
set(gca,'
XTick'
0:
2:
16)
7
XTickLabel'
{'
0架'
2架'
4架'
6架'
8架'
10架'
12架'
14架'
16架'
})
13
0.2460.2512
110.2110.20.210
90.1580.120.113probability70.16
0.05550.050.028140.01970.005533002468101214160架2架4架6架8架10架12架14架16架Numberofaircraftlaunchedontime
图10图11根据例5中的数据,用线性回归的方法画出95%的置信区间,见图11,程序如下:
x=[1:
[p,resid1]=polyfit(x,y,2);
[yhat,w]=polyconf(p,x,resid1,0.05);
plot(x,yhat,'
k-'
x,yhat-w,'
k--'
x,yhat+w,'
x,y,'
ks'
[x;
x],[yhat;
y],'
残差的进一步分析:
首先计算残差值,然后用函数normplot将其画在图12中,以观察残差是否服从正态分布。
程序为:
normplot(y-polyval(p,x))
whitebg('
white'
)%设置背景为白色
90
NormalProbabilityPlot800.98
700.95
0.9060
0.7550
400.50Probability300.25
200.10100.05
0.020155160165170175180185190195-1-0.500.5Data
图12图13由图12看出,这些残差点与直线非常接近,由此得出结论,这些残差值非常接近于正态分布,选择的模型是合适的。
8