MATLAB应用基础第四章.docx
- 文档编号:7844542
- 上传时间:2023-01-26
- 格式:DOCX
- 页数:39
- 大小:395.09KB
MATLAB应用基础第四章.docx
《MATLAB应用基础第四章.docx》由会员分享,可在线阅读,更多相关《MATLAB应用基础第四章.docx(39页珍藏版)》请在冰豆网上搜索。
MATLAB应用基础第四章
第4章MATLAB的数值计算
4.1求解线性代数方程组
已知线性代数方程组
Ax=b
可用下面三种方法求解
1)直接求逆法求解
x=inv(A)*b
此法当A不可逆时失效
2)左除法求解
x=A\b
左除法的基础是高斯消元法,由消元法对系数矩阵A进行LU分解。
进而得到方程组的解。
该法运算量少,运算速度快,而且数值稳定性好,解的精度高。
当方程个数大于未知变量个数时,该法可得到问题的最小二乘解。
当方程个数小于未知变量个数时,该法可求得有多个0元素的解。
3)使用伪逆函数求最小范数解:
x=pinv(A)*b
例如:
%3个未知量4个方程
A=[123;456;780;258]
b=[366,804,351,514]'
%计算最小二乘解
x=A\b
%该馀向量具有最小范数
res=A*x-b
%生成4个未知量的3个方程
A=A'
b=b(1:
3)
%具有最多0元素
x=A\b
%计算最小范数解
xn=pinv(A)*b
运行结果如下:
A=
123
456
780
258
b=
366
804
351
514
x=
247.9818
-173.1091
114.9273
res=
-119.4545
11.9455
0
35.8364
A=
1472
2585
3608
b=
366
804
351
x=
0
-165.9000
99.0000
168.3000
xn=
30.8182
-168.9818
99.0000
159.0545
4.2多项式计算
1、多项式的表示方法
多项式使用按降幂排列的多项式系数所构成的行向量描述。
即C1Xn+C2Xn–1+…+CnX+Cn+1可描述为:
[C1,C2,…,Cn,Cn+1]
2、多项式相乘
conv(a,b)
其中:
a,b为多项式系数向量
3、多项式相除
[q,r]=deconv(a,b)
其中:
q为商,r为余项
4、多项式相加
两个相同规模的多项式可以相加:
a+b
注意:
若规模不同,则要把低阶多项式补足若干个"0"。
5、多项式微分
1)求多项式p的导数
polyder(p)
2)求多项式之积的导数
polyder(a,b)
3)求多项式之比
的导数
[q,d]=polyder(a,b)
6、多项式求值
polyval(p,x)
功能:
求多项式p对应于x的值
7、多项式求根
roots(p)
8、求一组根所对应的多项式
若已知一多项式的根为行向量r时,则
poly(r)
为r所对应的多项式。
例如:
a=[1234];
b=[14916];
c=conv(a,b)
d=a+b
e=c+[000d]
[q,r]=deconv(e,a)
h=polyder(c)
g=polyder(a,b)
[q,d]=polyder(a,b)
p=[1,-12,0,25,116];
v=polyval(p,2.5)
r=roots(p)
运行结果如下:
c=
162050758464
d=
261220
e=
162052819684
q=
14918
r=
00002612
h=
6308015015084
g=
6308015015084
q=
212423212
d=
1834104209288256
v=
30.0625
r=
11.7473
2.7028
-1.2251+1.4672i
-1.2251-1.4672i
4.3数值逼近
1、曲线拟合
若已知离散数据向量x,y,则
polyfit(x,y,N)
将采用最小二乘法构造一个N阶多项式
例如:
x=[0.1.2.3.4.5.6.7.8.91];
y=[-.447,1.978,3.28,6.16,7.08,7.34,
7.66,9.56,9.48,9.3,11.2];
p=polyfit(x,y,2)
x1=0:
.01:
1;
z=polyval(p,x1);
p10=polyfit(x,y,10)
z10=polyval(p10,x1);
plot(x,y,'bo-',x1,z,'k:
',x1,z10,'r')
legend('原始数据','2阶曲线拟和','10阶曲线拟和')
运行结果如下:
p=
-9.810820.1293-0.0317
p10=
1.0e+006*
Columns1through7
-0.46442.2965-4.87735.8233-4.29482.0211-0.6032
Columns8through11
0.1090-0.01060.0004-0.0000
2、一元插值
若已知离散数据向量x,y,则
yi=interp1(x,y,xi,'method')
将计算在xi处的函数值yi。
其中的参数method可为:
nearest最邻近插值
linear线性插值(缺省)
spline三次样条插值
cubic三次多项式插值(要求x等距)xi为向量时,yi也为向量。
例如:
x=0:
8;
y=sin(x);
x1=0:
.1:
8;
y1=interp1(x,y,x1,'nearest');
y2=interp1(x,y,x1,'linear');
y3=interp1(x,y,x1,'spline');
y4=interp1(x,y,x1,'cubic');
plot(x,y,'o-',x1,y1,'b--',x1,y2,'k:
',x1,y3,'r-',x1,y4,'g-')
legend('原始数据','nearest','linear','spline','cubic')
运行结果如下:
3、二元插值
若已知离散数据x,y,z,则
①zi=interp2(x,y,z,xi,yi,'method')
功能:
计算在xi,yi处的zi值。
其中:
method可为:
linear线性插值(x,y单增)
bilinear双线性插值(x,y单增)
cubic三次插值(要求x,y单
增等距)
bicubic双三次插值(要求x,y
单增等距)
nearest最近邻域插值(要求x,
y单增)
②zi=griddata(x,y,z,xi,yi)
功能:
griddata使用距离倒数的方法,计算非等距节点x,y,z在xi,yi处的zi值。
若xi,yi为向量,则zi将为矩阵。
(xi为行向量,yi为列向量)。
xi和yi也可为由meshgrid函数生成矩阵。
例如:
%在500米间距正方形网格系统上测得的海底深度(米)
%绘制海底深度图
x=0:
0.5:
4;%0.5km为单位
y=0:
0.5:
6;%0.5km为单位
z=[1009910099100999999100
100999999100991009999
9999989810099100100100
1009897979910010010099
1011009898100102103100100
102103101100102106104101100
9910210010010310810610199
9799100100102105103101100
10010210310110210310210099
1001021031021011011009999
1001001011011001001009999
10010010010010099999999
10010010099991009910099];
mesh(x,y,z)
xlabel('X-axis(km)')
ylabel('Y-axis(km)')
zlabel('海底深度(m)')
title('海底深度图')
pause
xi=linspace(0,4,30);
yi=linspace(0,6,40);
[xxi,yyi]=meshgrid(xi,yi);
zzi=interp2(x,y,z,xxi,yyi,'linear');
mesh(xxi,yyi,zzi)
title('海底深度图(linear)')
holdon
[xx,yy]=meshgrid(x,y);
plot3(xx,yy,z+0.1,'ob')
holdoff
pause
zzi=interp2(x,y,z,xxi,yyi,'bilinear');
mesh(xxi,yyi,zzi)
title('海底深度图(bilinear)')
holdon
[xx,yy]=meshgrid(x,y);
plot3(xx,yy,z+0.1,'ob')
holdoff
pause
zzi=interp2(x,y,z,xxi,yyi,'cubic');
mesh(xxi,yyi,zzi)
title('海底深度图(cubic)')
holdon
[xx,yy]=meshgrid(x,y);
plot3(xx,yy,z+0.1,'ob')
holdoff
pause
zzi=interp2(x,y,z,xxi,yyi,'bicubic');
mesh(xxi,yyi,zzi)
title('海底深度图(bicubic)')
holdon
[xx,yy]=meshgrid(x,y);
plot3(xx,yy,z+0.1,'ob')
holdoff
pause
zzi=interp2(x,y,z,xxi,yyi,'nearest');
mesh(xxi,yyi,zzi)
title('海底深度图(nearest)')
holdon
[xx,yy]=meshgrid(x,y);
plot3(xx,yy,z+0.1,'ob')
holdoff
pause
zzi=interp2(x,y,z,xxi,yyi,'cubic');
pcolor(xxi,yyi,zzi)
shadinginterp
holdon
contour(xxi,yyi,zzi,15,'k')
colormap(cool)
colorbar('vert')
holdoff
title('海底深度等值线图')
运行结果如下:
4.4常微分方程的数值解
由于任一高阶常微分方程均可化为与之等价的一阶常微分方程组(即"状态方程"),所以MATLAB中总是求解一阶常微分方程组,而且需要用m文件定义方程组。
例如:
若常微分方程为:
设
则原方程可表为
故可定义xp.m文件描述为:
functionxdot=xp(t,y)
globalMU%0 xdot=zeros(2,1); xdot (1)=y (2); xdot (2)=MU*(1–y (1)^2)*y (2)–y (1); 该函数对已知的变量返回列向量 注意: 在函数xp中,自变量必须给出两个,头一个必须是"时间变量",第二个必须是方程组向量("状态变量")。 假定: t的变化区间为[0,50], y0的初值为(0,0.25) 则可用下述命令求解: [t,y]=0de23('xp',[0,20],[0;0.25]) 或者 [t,y]=0de45('xp',[0,20],[0;0.25]) 其中: [t,y]=0de23('Fun',t,y0) 功能是: 在指定区间[t0,tf]上,使用2阶和3阶龙格-库塔法求解满足初始条件列向量Y0的、由函数Fun.m所描述的常微分方程组,返回时间列向量t和矩阵y,t的每个元素对应y的一行。 [t,y]=0de23('Fun',t,y0,to,trac) 功能同上,可指定求解的精度to(缺省值为相对误差1e-3,绝对误差1e-6),当trac非0时,将输出每步中间的结果(缺省为0)。 [t,y]=0de45('Fun',t,y0) [t,y]=0de45('Fun',t,y0,to,trac) 功能: 使用4阶和5阶龙格-库塔法求解 另外,ode23s、ode23t、ode23tb函数,可给出解的动态轨迹;ode113用可变阶方法解非Stiff方程;ode15s用可变阶方法解Stiff方程。 例如: t=[0,50]; y0=[0,0.25]; [t,y]=ode45('xp',t,y0); plot(t,y(: 1),': b',t,y(: 2),'-r') legend('速度','位移') pause clf; axis([0,50,-4,4]); holdon ode23s('xp',t,y0) legend('速度','位移') 运行结果如下: 4.5数值积分 MATLAB中有多个求定积分值的函数: 1.quad('fun',a,b) 功能: 使用自适应辛普生法求函数fun在[a,b]上的定积分值。 2.quad('fun',a,b,to,trac) 功能: 使用自适应辛普生法求函数fun在[a,b]上的定积分值,可指定求解的精度to(缺省值为相对误差1e-3),当trac非0时,将画出被积函数的点迹(缺省为0)。 3.quad8('fun',a,b) 功能: 使用自适应递归的牛顿-考梯斯法(Newton-Coties)计算函数fun在[a,b]上的定积分值。 其中: fun必须是一个fun.m文件。 4.trapz(x,y) 功能: 使用梯形公式计算定积分值。 被积函数用向量x和y给出。 4.6方程求根 1、求一元函数的零点 x=fzero('fun',x0) 功能: 将求得自x0开始的函数fun的一个零点。 其中fun通常是文件fun.m的名字,也可以是一个字符串。 例如: x=fzero('sin',3) 的执行结果为: Zerofoundintheinterval: [2.8303,3.1697]. x= 3.1416 2、解非线性方程组 x=fsolve('fun',x0) 或 x=fsolve('fun',x0,options) 功能: 将用最小二乘法求得自x0开始的非线性方程组fun的一个零点。 其中fun是文件fun.m的名字 其中: options为一参数向量,通过该向量可设定控制参数。 运行helpfoptions命令的部分结果如下: OPTIONS (1)-Displayparameter(Default: 0). 1displayssomeresults. OPTIONS (2)-TerminationtoleranceforX. (Default: 1e-4). OPTIONS(3)-TerminationtoleranceonF. (Default: 1e-4). OPTIONS(4)-Terminationcriteriononconstraint violation.(Default: 1e-6) OPTIONS(5)-Algorithm: Strategy: Notalwaysused. OPTIONS(6)-Algorithm: Optimizer: Notalwaysused. OPTIONS(7)-Algorithm: LineSearchAlgorithm. (Default0) OPTIONS(8)-Functionvalue.(Lambdaingoal attainment.) OPTIONS(9)-Setto1ifyouwanttocheck user-suppliedgradients OPTIONS(10)-NumberofFunctionandConstraint Evaluations. OPTIONS(11)-NumberofFunctionGradient Evaluations. OPTIONS(12)-NumberofConstraintEvaluations. OPTIONS(13)-Numberofequalityconstraints. OPTIONS(14)-Maximumnumberoffunction evaluations.(Defaultis100*numberof variables) OPTIONS(15)-Usedingoalattainmentforspecial objectives. OPTIONS(16)-Minimumchangeinvariablesfor finitedifferencegradients. OPTIONS(17)-Maximumchangeinvariablesfor finitedifferencegradients. OPTIONS(18)-Steplength.(Default1orless). 例如: 若方程组为: 可用以下xyz.m文件描述它: functionq=xyz(p) x=p (1); y=p (2); z=p(3); q=zreos(3,1); q (1)=sin(x)+y^2+log(z)-7; q (2)=3*x+2^y-z^3+1; q(3)=x+y+z-5; 下面指定初始点X0,并求解: x0=[1,1,1]; x=fsolve('xyz',x0) 运行结果如下: x= 0.59952.39592.0051 4.7求解无约束的最优化问题 无约束的最优化问题一般形式为: minf(x) 1.求一元函数的极小点(区间为[x1,x2]) x=fmin('fun',x1,x2) 2.用单纯形法求多元函数在x0点附近的极小点 x=fmins('fun',x0) 3、用BFGS拟牛顿法由多元函数在x0点附近的极小点 x=fminu('fun',x0) 4、用最小二乘法求多元函数在x0点附近的极小点 x=leastsq('fun',xo) 例如: %BANDEMOBananafunctionminimizationdemonstration. % %Thisscript-filedemonstratestheminimizationofthe %thebananafunction: % %f(x)=100*(x (2)-x (1)^2)^2+(1-x (1))^2 % %Itiscalledthebananafunctionbecauseofthewaythe %curvaturebendsaroundtheorigin(alsocalledRosenbrock's %function).Itisnotoriousinoptimizationexamplesbecause %oftheslowconvergencewithwhichmostmethodsexhibit %whentryingtosolvethisproblem. % %ThisfunctionhasauniqueminimumatthepointX=[1,1]f(x)=0. %Wedemonstratehereanumberoftechniquesforitsminimization %startingatthepointX=[-1.9,2]. % echooff clc helpbandemo disp('Strikeanykeyforameshplotofthebananafunction') xx=[-2: 0.125: 2]'; yy=[-1: 0.125: 3]'; [x,y]=meshgrid(xx',yy'); meshd=100.*(y-x.*x).^2+(1-x).^2; pause mesh(meshd) disp('') disp('Strikeanykeyforacontourplotofthebananafunction') pause clf conts=exp(3: 20); contour(xx,yy,meshd,conts,'w--') xlabel('x1') ylabel('x2') title('MinimizationoftheBananafunction') drawnow;%Drawscurrentgraphnow holdon plot(-1.9,2,'o') text(-1.9,2,'StartPoint') plot(1,1,'o') text(1,1,'Solution') disp('Pleasewait-compilingoptimizationroutines') %test_longisavariableusedforautotestingofthisroutine if~exist('test_long')test_long=0;end ifexist('method')~=1method=7;end if~length(method)method=7;end l=2; while1 if~test_long clc disp('') disp('Chooseanyofthefollowingmethodstominimizethebananafunction') disp('') disp('UNCONSTRAINED: 1)Broyden-Fletcher-Golfarb-Shanno') disp('2)Davidon-Fletcher-Powell') disp('3)SteepestDescent') disp('4)SimplexSearch') disp('LEASTSQUARES: 5)Gauss-Newton') disp('6)Levenberg-Marqua
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 应用 基础 第四