第六章MATLAB 数值计算.docx
- 文档编号:7759224
- 上传时间:2023-01-26
- 格式:DOCX
- 页数:26
- 大小:259.80KB
第六章MATLAB 数值计算.docx
《第六章MATLAB 数值计算.docx》由会员分享,可在线阅读,更多相关《第六章MATLAB 数值计算.docx(26页珍藏版)》请在冰豆网上搜索。
第六章MATLAB数值计算
第六章MATLAB数值计算
6-1多项式的运算
6-1-1多项式的生成和表达
1.多项式的表达
在MATLAB环境下多项式是用向量的形式表达的.向量最右边的元素表示多项式的0阶,向左数依次表示多项式的第1阶、第2阶、第3阶…。
例如多项式
表示为:
[50321]。
2.多项式的生成
语法:
P=ploy(MA)
说明:
1.若MA为方阵,则生成的多项式P为方阵MA的特征多项式。
2.若MA为向量,则向量和多项式满足这样一种关系:
,生成的多项式为:
3.直接输入的方式生成多项式。
例6-1
利用方阵M=[567;891;111213]生成一个多项式(为方阵M的特征多项式)。
程序设计:
>>clear
M=[567;891;111213];
P=poly(M);%产生多项式的向量表达式
Px=poly2str(P,'x');%生成常见的多项式表示形式
P,Px
运行结果:
P=
1.0000-27.000090.000054.0000
Px=
x^3-27x^2+90x+54
例6-2
利用向量A=[2345]生成一个多项式。
程序设计:
>>clear
A=[2345];
P=poly(A);%产生多项式的向量表达式
Px=poly2str(P,'x');%生成常见的多项式表示形式
P,Px
运行结果:
P=
1-1471-154120
Px=
x^4-14x^3+71x^2-154x+120
6-1-2多项式的乘除
语法:
A.c=conv(a,b)
B.[q,r]=decony(c,a)
说明:
1.a、b和c分别是多项式的向量表示形式。
A表示两个多项式的乘积运算,B表示两个多项式的除法运算。
2.q表示除运算的商,r表示除运算的余数。
例6-3
求多项式
和
的乘积M(x)。
程序设计:
>>clear
a=[150];%第一个多项式F(x)
b=[21];%第二个多项式G(x)
c=conv(a,b);%求两个多项式的乘积
Mx=poly2str(c,'x');%用常用的方式表示多项式的积
c,Mx
%end
运行结果:
c=
21150
Mx=
2x^3+11x^2+5x
例6-4
求多项式
和
的除运算D(x)。
程序设计:
>>clear
c=[150];%第一个多项式F(x)
a=[21];%第二个多项式G(x)
[q,r]=deconv(c,a);%求F(x)/G(x)
Dx=poly2str(q,'x');%用常用的方式表示多项式的积
q,r,Dx
%end
运行结果:
q=
0.50002.2500
r=
00-2.2500
Dx=
0.5x+2.25
程序说明:
1.在运行结果中变量q是F(x)除以G(x)的商,而r则是除不尽的余数。
2.运行结果变量Dx表示的商没有加上余数。
6-1-3多项式的求导
语法:
Dp=polyder(p)
说明:
p为向量表示的多项式。
例6-5
求多项式
和
的一阶导数。
我们容易知道以上两个方程的导数手工验算结果为:
F'(x)=2x+5和G'(x)=2
我们看MATLAB的计算结果。
程序设计:
>>clear
f=[150];%第一个多项式F(x)
g=[21];%第二个多项式G(x)
Df=polyder(f);%求F(x)的导数
Dg=polyder(g);%求G(x)的导数
Dfx=poly2str(Df,'x');
Dgx=poly2str(Dg,'x');
Df,Dg,Dfx,Dgx
%end
运行结果:
Df=
25
Dg=
2
Dfx=
2x+5
Dgx=
2
6-1-4多项式的求根
语法:
A.r=roots(p)
说明:
另外还有一种通过先求多项式的伴随矩阵,然后再求特征值的方法也可以求得多项式的根。
例6-6
求多项式
和
的根。
我们容易知道方程的根为:
F(x)为:
x1=0;x2=-5
G(x)为:
x2=-1/2
我们看MATLAB的计算结果。
程序设计:
>>clear
f=[150];%第一个多项式F(x)
g=[21];%第二个多项式G(x)
rf=roots(f);%求F(x)的根
rg=roots(g);%求G(x)的根
rf,rg,
%end
运行结果:
rf=
0
-5
rg=
-0.5000
程序说明:
从运行结果我们可以看到F(x)的根为0和-5,G(x)的根为-0.50,这和我们手工验算的结果完全一致。
例6-7
求多项式
和
的根。
程序设计:
>>clear
f=[15-3];%第一个多项式F(x)
g=[2101];%第二个多项式G(x)
rf=roots(f);%求F(x)的根
rg=roots(g);%求G(x)的根
rf,rg,
%end
运行结果:
rf=
-5.5414
0.5414
rg=
-1.0000
0.2500+0.6614i
0.2500-0.6614i
程序说明:
1.我们可以看到例6-6和例6-7求得的根和我们手工计算的结果是一致的,其中例6-7多项式的根出现了虚根。
2.我们用求伴随矩阵的方法对例6-7再做一次求根运算:
>>clear
f=[15-3];%第一个多项式F(x)
g=[2101];%第二个多项式G(x)
cf=compan(f);%求F(x)的伴随矩阵
cg=compan(g);%求G(x)的伴随矩阵
cf,cg,
%end
运行结果:
cf=
-53
10
cg=
-0.50000-0.5000
1.000000
01.00000
ccf=eig(cf)%求特征根
ccg=eig(cg)%求特征根
运行结果:
ccf=
-5.5414
0.5414
ccg=
-1.0000
0.2500+0.6614i
0.2500-0.6614i
我们看到两种方法的结果是一致的。
6-2数据分析
6-1-1极值、均值、标准差和中位值的计算
语法:
A.Pmax=max(X)
B.Pmin=min(X)
C.Pmean=mean(X)
D.Pstd=std(X)
E.Pmedian=median(X)
说明:
1.max(X)、min(X)、mean(X)、std(X)、median(X)分别是用来求数组或矩阵的最大值、最小值、均值、标准差、中位值。
注意std(X)的定义是
,std(X,1)的定义是
。
2.这里X可以是数组也可以是矩阵。
如果是数组则对整个数组进行计算,如果是矩阵则是分别对矩阵的每个列向量进行运算。
例6-8
对一组随机数组进行均值、方差和中位值的计算。
程序设计:
>>clear
x=randn(1,10);%生成一个10元素的数组
Pmax=max(x);%计算数组x的最大值
Pmin=min(x);%计算数组x的最小值
Pmean=mean(x);%计算数组x的平均值
Pstd=std(x);%计算数组x的标准差
Psqu=Pstd^2;%计算数组x的方差
Pmed=median(x);%计算数组x的中位值
answers=[PmaxPminPmeanPstdPsquPmed];
x,answers
运行结果:
x=
Columns1through6
-0.4326-1.66560.12530.2877-1.14651.1909
Columns7through10
1.1892-0.03760.32730.1746
answers=
1.1909-1.66560.00130.90340.81620.1500
程序说明:
1.注意方差为标准差的平方。
2.需要注意的是MATLAB的变量定义是区分大小写的。
例6-9
对一个随机矩阵进行最大值、最小值、均值、标准差、中位值的计算。
程序设计:
>>clear
X=randn(6,3);%生成一个6×3的随机矩阵
Pmax=max(X);%计算矩阵X的最大值
Pmin=min(X);%计算矩阵X的最小值
Pmean=mean(X);%计算矩阵X的平均值
Pstd=std(X);%计算矩阵X的标准差
Pmed=median(X);%计算矩阵X的中位值
X,Pmax,Pmin,Pmean,Pstd,Pmed;
运行结果:
X=
-0.18671.06680.7143
0.72580.05931.6236
-0.5883-0.0956-0.6918
2.1832-0.83230.8580
-0.13640.29441.2540
0.1139-1.3362-1.5937
Pmax=
2.18321.06681.6236
Pmin=
-0.5883-1.3362-1.5937
Pmean=
0.3519-0.14060.3607
Pstd=
0.99620.84821.2404
6-2-2曲线拟合
语法:
A.P=polyfit(X,Y,N)
B.Yval=polyval(P,X)
说明:
1.polyfit(X,Y,N)根据输入数据X和Y生成一个N阶的拟合多项式。
2.polyval(P,X)根据数据X,用拟合多项式P生成拟合好的数据。
3.这里X和Y构成了一系列数据点的坐标。
例6-10
下列数据为一段时间的日气温平均实际数据,求该数据的拟合方程。
设这组气温数据为:
[181917182022222321212022232322]
程序设计:
>>clear
d=1:
15;%时间
t=[181917182022222321212022232322];%和时间相对应的日平均气温
p=polyfit(d,t,3);%生成拟合多项式(向量形式),阶次为3
px=poly2str(p,'x');%生成拟合多项式(字符形式)
pv=polyval(p,d);%生成拟合值
p,px
plot(d,t,':
*r',d,pv,'-k','LineWidth',2)
title('\fontsize{18}\bf数据拟合图像','Color','r')
xlabel('\fontsize{14}\rm时间d轴','Color','k')
ylabel('\fontsize{14}\rm气温t轴','Color','k')
gtext('\fontsize{16}\fontname{宋体}实际曲线','Color','r')
gtext('\fontsize{16}\fontname{宋体}拟合曲线','Color','k')
%end
运行结果:
p=
0.0010-0.05380.964416.4623
px=
0.x^3-0.x^2+0.96436x+16.4623
图形如6-1所示:
程序说明:
1.d表示时间第一天、第二天…,t表示每天的气温。
2.本例我们选取的拟合多项式的阶次是3。
3.p为向量形式的拟合多项式,px为字符形式的拟合多项式。
4.plot为绘图函数,它的两个参数分别为图像的横坐标和纵坐标。
4个参数则分别为两条曲线的横坐标和纵坐标。
5.图中的箭头可在绘图窗口中打开insert菜单,选择Arrow然后用鼠标从起始位置拖动到终止位置即可。
例6-11
根据上例中的数据建立的天气温度,拟合方程求5天后的天气气温是多少?
程序设计及运行结果:
>>clear
x=20;
T=0.*x^3-0.*x^2+0.96436*x+16.4623
T=
22.5995
程序设计及运行结果:
1.依照上例7-10天的时间序列d,5天后应该是20天,所以令x=20。
2.T即为上例中我们求得的拟合公式,结果求得5天后的日平均气温为22.5995。
6-3数值积分和数值微分
6-3-1微分和积分的数学表达式
对函数f(x)在区间[x1,x2]上求积分在数学上常常表示为
对函数F(x)求导数运算在数学上常常表示为:
6-3-2函数数值积分
命令函数:
1.连续被积函数
quadl('f',a,b,t,trace)自适应递推牛顿—科西(Newton—Cotes)法求积分
2.离散被积函数
trapz(X,Y)梯形法求数值积分
说明:
1.
1)参数'f'是被积函数表达式字符串或者函数文件名。
2)参数a、b定义函数积分的上下限。
3)参数t定义积分的精度(默认值1.0e-6);trace设置是否用图形展示积分过程,1为展示、0不展示。
2.
trapz积分中给出Y相对于X的积分值。
当Y是(m*n)矩阵时,积分对Y的列向量分别进行,得到一个(1*n)矩阵是Y的列向量对应于X的积分结果。
例6-12
设
,用牛顿—科西法求
。
这个积分可以用手工算出来是2,我们看MATLAB的计算结果。
程序设计及运行结果:
>>clear
s=quadl('sin(x)',0,pi)
s=
2.0000
例6-13
设
,求
。
程序设计及运行结果:
>>clear
y=quadl('(sin(x).^2+x.^3)./(x-2.*cos(3.*x))',0,2)
%end
ans=
11.1943
例6-14
用离散数据表示
,积分区间为
分别用离散积分方法对其进行积分。
程序设计及运行结果(如图6-3):
>>clear
X=0:
0.001:
pi;%离散化积分变量
Y=sin(X);%生成被积函数的离散序列
Z=trapz(X,Y);%梯形法求积分
Z
%end
Z=
2.0000
例6-15
设
,求
。
程序设计及运行结果:
>>clear
X=0:
0.001:
2;%离散化积分变量
Y=(sin(X).^2+X.^3)./(X-2.*cos(3.*X));
Z=trapz(X,Y);
Z
%end
Z=
10.9607
程序说明:
输入公式时一定要注意算符都是对数组进行运算的,所以是“.*”“.^”“./“。
6-3-3数值微分
命令函数:
D=diff(Y)
说明:
Y为一组离散序列,D为与Y同维的离散矩阵,则微分运算是对矩阵的每个列向量微分。
例6-16
设函数
和
进行微分运算,并画出计算结果的图像。
程序设计:
>>clear
dx=0.001;
x=-3:
dx:
3;%生成一个从-3到3间隔为0.001的序列
x1=-3:
dx:
3-dx;
y=sin(x);
z=2.*x.^3+3.*x.^2+x-1;
Y=diff(y);%对y进行微分运算
DY=diff(y)/dx;%对y进行求导运算
Z=diff(z);%对z进行微分运算
DZ=diff(z)/dx;%对z进行求导运算
subplot(2,2,1);%产生两行一列的绘图区间的第一个绘图区间
plot(Y,'-k','LineWidth',1)
title('\fontsize{12}\fontname{TimeNewRoman}\Deltay(x)','Color','k')
gridon
subplot(2,2,2);%产生两行一列的绘图区间的第二个绘图区间
plot(Z,'-k','LineWidth',1)
title('\fontsize{12}\fontname{TimeNewRoman}\Deltaz(x)','Color','k')
gridon
subplot(2,2,3);%产生两行一列的绘图区间的第三个绘图区间
plot(x1,DY,'-k','LineWidth',1)
title('\fontsize{12}\fontname{TimeNewRoman}y''(x)=cosx','Color','k')
xlabel('\fontsize{12}\fontname{TimeNewRoman}x','Color','k')
ylabel('\fontsize{12}\fontname{TimeNewRoman}dy/dx','Color','k')
axis([-33-11])
gridon
subplot(2,2,4);%产生两行一列的绘图区间的第四个绘图区间
plot(x1,DZ,'-k','LineWidth',1)
title('\fontsize{12}\fontname{TimeNewRoman}z''(x)=6x^2+6x+1','Color','k')
xlabel('\fontsize{12}\fontname{TimeNewRoman}x','Color','k')
ylabel('\fontsize{12}\fontname{TimeNewRoman}dz/dx','Color','k')
axis([-33-2080])
gridon
%end
运行结果(如图6-2):
程序说明:
1.输入公式时一定要注意算符都是对数组进行运算的,所以是“.*”“.^”“./”。
2.Y=diff(y)和Z=diff(z)和依此画出的图像,纵坐标是两函数的微分,横坐标是数据序列的值并不是实际图像的横坐标。
由于纵坐标是
,所以纵坐标在数值上是导数的
倍。
3.导数图的纵坐标是
,n的取值范围应该是
。
因此
取值是
到
,而微分相应取值是
到
,所以两序列的数目相差1,所以画图时不能用plot(x,diff(y)/dx)。
6-4一般非线性方程组的数值解
命令函数:
x=fsolve('fun',x0,options)
说明:
1.'fun'是我们要求解的方程组,可以直接输入,不过由于方程组比较复杂,一般都先生成一个m文件,然后再进行调用。
2.x0是给出的这个方程组的初值解,我们可以随意地给出。
3.options是命令函数fsolve的参数设置项。
因为fsolve求解的过程本身就是一个优化的过程,而options则是设置优化过程的参数的。
这里对于一般的非线性方程组求解我们常设置为optiomset('Display','off'),意思为不显示优化过程。
例:
6-17
求方程组
的数值解。
这个方程组的解经过手工演算我们知道解为:
x=4/3,y=-2/3,z=-2/3。
程序设计与运行结果:
1.建立函数文件xyz,并保存:
functionfun=xyz(X)
x=X
(1);y=X
(2);z=X(3);
fun=zeros(3,1);
fun
(1)=x+y+z;
fun
(2)=x-y+3*z;
fun(3)=3*x+y-z-4;
2.在MATLAB的工作区里求解:
>>X0=[111];
X=fsolve(@xyz,X0,optimset('Display','off'))
X=
1.3333-0.6667-0.6667
程序说明:
1.在建立一个函数文件时一定要以function作为文件头说明。
2.把建立的函数文件保存到MATLAB默认的文件夹里,如work文件夹。
3.zeros(3,1)用来产生一个3行1列的矩阵,矩阵的所有元素都为0。
4.我们可以看到计算的结果和我们手工计算的结果非常近似,差别是计算精度引起的,fsolve的默认求解精度是0.001。
5.事实上对方程组的求解过程是一个优化过程,optimset('Display','off')设置不显示求解过程。
例:
6-18
求非线性方程组
的数值解。
这个方程组手工计算是难以完成的。
程序设计与运行结果:
1.建立函数文件fun,并保存:
functionQ=fun(T)
x=T
(1);y=T
(2);z=T(3);
Q=zeros(3,1);
Q
(1)=x+y-z;
Q
(2)=cos(x)+y^3+z-5;
Q(3)=4*x+2^y+log(z)-1;
2.在MATLAB的工作区里求解:
>>clear
X=fsolve(@fun,[111],optimset('Display','off'))
X=
-0.43921.45481.0156
程序说明:
fsolve的默认求解精度是0.001,方程成立的精度为0.001。
6-5微分方程的数值解
6-5-1微分方程的基本形式
的导数等于
,即:
6-5-2一阶常微分方程的求解
语法:
A.[T,Y]=ode45(odefun,tspan,y0)
B.[T,Y]=ode23(odefun,tspan,y0)
说明:
1.ode45()和ode23()是两个求解微分方程的命令函数,ode45()应用最广,在容许误差大的情况下ode23()的效率要比ode45()高一些。
2.T和Y分别表示微分方程的解的两个变量的数值序列。
3.odefun是待解微分方程表达式的函数文件名。
4.tspan表示运算的起至时间,是行向量,例如[010]表示运算时间从0开始到10结束。
5.y0初始状态,用列向量表示,其元素分别表示微分方程组各个表达式的初始状态值。
例:
6-19
求解微分方程
并画出解在区间[02*pi]上的图像。
程序设计与运行结果:
1.在程序编辑器中建立待解微分方程表达式的函数文件funx.m,并存入work工作目录:
functionY=funx(x,y)
Y=cos(x);
2.在MATLAB的命令窗口求解,结果如图6-3所示:
>>clear
[x,y]=ode45(@funx,[02*pi],0);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第六章 MATLAB 数值计算 第六 数值 计算