第06章MATLAB数值计算例题源程序汇总.docx
- 文档编号:27627218
- 上传时间:2023-07-03
- 格式:DOCX
- 页数:17
- 大小:30.05KB
第06章MATLAB数值计算例题源程序汇总.docx
《第06章MATLAB数值计算例题源程序汇总.docx》由会员分享,可在线阅读,更多相关《第06章MATLAB数值计算例题源程序汇总.docx(17页珍藏版)》请在冰豆网上搜索。
第06章MATLAB数值计算例题源程序汇总
第6章MATLAB数值计算
例6.1求矩阵A的每行及每列的最大和最小元素,并求整个矩阵的最大和最小元素。
A=[13,-56,78;25,63,-235;78,25,563;1,0,-1];
max%求每行最大元素
min%求每行最小元素
max%求每列最大元素
min%求每列最小元素
max
也可使用命令:
max >> min 也可使用命令: min >> 例6.2求矩阵A的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; prod 也可以使用命令prod >> 例6.3求向量X=<1! 2! 3! …,10! >。 X=cumprod<1: 10> 例6.4对二维矩阵x,从不同维方向求出其标准方差。 x=[4,5,6;1,4,8]%产生一个二维矩阵x y1=std y2=std y3=std y4=std 例6.5生成满足正态分布的10000×5随机矩阵,然后求各列元素的均值和标准方差,再求这5列随机数据的相关系数矩阵。 X=randn<10000,5>; M=mean D=std R=corrcoef 例6.6对下列矩阵做各种排序。 A=[1,-8,5;4,12,6;13,7,-13]; sort%对A的每列按升序排序 -sort<-A,2>%对A的每行按降序排序 [X,I]=sort%对A按列排序,并将每个元素所在行号送矩阵I 例6.7给出概率积分 的数据表如表6.1所示,用不同的插值方法计算f<0.472>。 表6.1概率积分数据表 x 0.46 0.47 0.48 0.49 f 0.4846555 0.4937542 0.5027498 0.5116683 x=0.46: 0.01: 0.49;%给出x,f f=[0.4846555,0.4937542,0.5027498,0.5116683]; formatlong interp1 interp1 interp1 interp1 formatshort 例6.8某检测参数f随时间t的采样结果如表6.2,用数据插值法计算t=2,7,12,17,22,17,32,37,42,47,52,57时的f值。 表6.2检测参数f随时间t的采样结果 t 0 5 10 15 20 25 30 f 3.1025 2.256 879.5 1835.9 2968.8 4136.2 5237.9 t 35 40 45 50 55 60 65 f 6152.7 6725.3 6848.3 6403.5 6824.7 7328.5 7857.6 T=0: 5: 65; X=2: 5: 57; F=[3.2015,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,... 6725.3,6848.3,6403.5,6824.7,7328.5,7857.6]; F1=interp1 F1=interp1 F1=interp1 F1=interp1 例6.9设z=x2+y2,对z函数在[0,1]×[0,2]区域内进行插值。 x=0: 0.1: 1;y=0: 0.2: 2; [X,Y]=meshgrid Z=X.^2+Y.^2;%求对应的函数值 interp2 interp2 interp2 %下一命令在<0.5,0.4>,<0.6,0.4>,<0.5,0.5>和<0.6,0.5>各点插值 interp2 例6.10某实验对一根长10米的钢轨进行热源的温度传播测试。 用x表示测量点<米>,用h表示测量时间<秒>,用T表示测得各点的温度<℃>,测量结果如表6.2所示。 表6.3钢轨各点温度测量值 Tx h 02.557.510 0 9514000 30 884832126 60 6764544841 试用3次多项式插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。 x=0: 2.5: 10; h=[0: 30: 60]'; T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41]; xi=[0: 0.5: 10]; hi=[0: 10: 60]'; temps=interp2 mesh 例6.11用一个3次多项式在区间[0,2π]内逼近函数 。 X=linspace<0,2*pi,50>; Y=sin P=polyfit X=linspace<0,2*pi,20>; Y=sin Y1=polyval plot o',X,Y1,'-*'> 例6.12设 <1>求f <2>求f f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g]; f+g1%求f f-g1%求f conv [Q,r]=deconv 例6.13求有理分式的导数。 P=[3,5,0,-8,1,-5]; Q=[10,5,0,0,6,0,0,7,-1,0,-100]; [p,q]=polyder 例6.14已知多项式x4+8x3-10,分别取x=1.2和一个2×3矩阵为自变量计算该多项式的值。 A=[1,8,0,0,-10];%4次多项式系数 x=1.2;%取自变量为一数值 x=[-1,1.2,-1.4;2,-1.8,1.6]%给出一个矩阵x y2=polyval%分别计算矩阵x中各元素为自变量的多项式之值 例6.15仍以多项式x4+8x3-10为例,取一个2×2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。 A=[1,8,0,0,-10];%多项式系数 x=[-1,1.2;2,-1.8]%给出一个矩阵x y1=polyval%计算代数多项式的值 y2=polyvalm%计算矩阵多项式的值 例6.16求多项式x4+8x3-10的根。 A=[1,8,0,0,-10]; 例6.17已知 <1>计算f <2>由方程f P=[3,0,4,-5,-7.2,5]; X=roots %求方程f G=poly 例6.18设x由[0,2π]间均匀分布的10个点组成,求 的1~3阶差分。 X=linspace<0,2*pi,10>; Y=sin DY=diff D2Y=diff D3Y=diff 例6.19设 用不同的方法求函数f f=inline<'sqrt g=inline<'<3*x.^2+4*x-1>./sqrt x=-3: 0.01: 3; p=polyfit dp=polyder ;%对拟合多项式p求导数dp dpx=polyval dx=diff gx=g plot 例6.20用两种不同的方法求: 先建立一个函数文件ex.m: functionex=ex ex=exp<-x.^2>; 然后在MATLAB命令窗口,输入命令: formatlong I=quad<'ex',0,1>%注意函数名应加字符引号 I=quadl<'ex',0,1> 例6.21用trapz函数计算: X=0: 0.01: 1; Y=exp<-X.^2>; trapz 例6.22计算二重定积分 <1>建立一个函数文件fxy.m: functionf=fxy globalki; ki=ki+1;%ki用于统计被积函数的调用次数 f=exp<-x.^2/2>.*sin <2>调用dblquad函数求解。 globalki;ki=0; I=dblquad<'fxy',-2,2,-1,1> ki 例6.23给定数学函数 x 取N=128,试对t从0~1秒采样,用FFT作快速傅立叶变换,绘制相应的振幅-频率图。 N=128;%采样点数 T=1;%采样时间终点 t=linspace<0,T,N>;%给出N个采样时间ti N> x=12*sin<2*pi*10*t+pi/4>+5*cos<2*pi*40*t>;%求各采样点样本值x dt=t<2>-t<1>;%采样周期 f=1/dt;%采样频率 X=fft F=X<1: N/2+1>;%F N/2+1> f=f*<0: N/2>/N;%使频率轴f从零开始 plot xlabel<'Frequency'>; ylabel<'|F 例6.24用直接解法求解下列线性方程组。 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b 例6.25用LU分解求解例6.24中的线性方程组。 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [L,U]=lu; x=U\ [L,U,P]=lu; x=U\ 例6.26用QR分解求解例6.24中的线性方程组。 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [Q,R]=qr; x=R\ [Q,R,E]=qr; x=E* 例6.27用Cholesky分解求解例6.24中的线性方程组。 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; ? ? ? Errorusing==>chol Matrixmustbepositivedefinite Jacobi迭代法的MATLAB函数文件Jacobi.m如下: ifnargin==3 eps=1.0e-6; elseifnargin<3 error return end D=diag L=-tril;%求A的下三角阵 U=-triu;%求A的上三角阵 B=D\ f=D\b; y=B*x0+f; n=1;%迭代次数 whilenorm x0=y; y=B*x0+f; n=n+1; end 例6.28用Jacobi迭代法求解下列线性方程组。 设迭代初值为0,迭代精度为10-6。 在命令中调用函数文件Jacobi.m,命令如下: A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]'; Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下: ifnargin==3 eps=1.0e-6; elseifnargin<3 error return end D=diag L=-tril;%求A的下三角阵 U=-triu;%求A的上三角阵 G= f= y=G*x0+f; n=1;%迭代次数 whilenorm x0=y; y=G*x0+f; n=n+1; end 例6.29用Gauss-Serdel迭代法求解下列线性方程组。 设迭代初值为0,迭代精度为10-6。 在命令中调用函数文件gauseidel.m,命令如下: A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]'; 例6.30分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组,看是否收敛。 a=[1,2,-2;1,1,1;2,2,1]; b=[9;7;6]; 有了上面这些讨论,下面设计一个求解线性方程组的函数文件line_solution.m。 [m,n]=size; y=[]; ifnorm>0%非齐次方程组 ifrank==rank<[A,b]> ifrank==n%有惟一解 disp<'原方程组有惟一解x'>; x=A\b; else%方程组有无穷多个解,基础解系 disp<'原方程组有无穷个解,特解为x,其齐次方程组的基础解系为y'>; x=A\b; y=null; end else disp<'方程组无解'>;%方程组无解 x=[]; end else%齐次方程组 disp<'原方程组有零解x'>; x=zeros ifrank disp<'方程组有无穷个解,基础解系为y'>;%非0解 y=null; end end 例6.31求解方程组。 A=[1,-2,3,-1;3,-1,5,-3;2,1,2,-2]; b=[1;2;3]; 例6.32求方程组的通解。 formatrat%指定有理式格式输出 A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8]; b=[1,4,0]'; [x,y]=line_solution; x,y formatshort%恢复默认的短格式输出 例6.33求 在x0=-5和x0=1作为迭代初值时的零点。 先建立函数文件fz.m: functionf=fz f=x-1/x+5; 然后调用fzero函数求根。 : fzero<'fz',-5>%以-5作为迭代初值 fzero<'fz',1>%以1作为迭代初值 例6.34求下列方程组在<1,1,1>附近的解并对结果进行验证。 首先建立函数文件myfun.m。 functionF=myfun x=X<1>; y=X<2>; z=X<3>; F<1>=sin F<2>=x+y+z; F<3>=x*y*z; 在给定的初值x0=1,y0=1,z0=1下,调用fsolve函数求方程的根。 X=fsolve<'myfun',[1,1,1],optimset<'Display','off'>> 例6.35求圆和直线的两个交点。 圆: 直线: 先建立方程组函数文件fxyz.m: functionF=fxyz x=X<1>; y=X<2>; z=X<3>; F<1>=x^2+y^2+z^2-9; F<2>=3*x+5*y+6*z; F<3>=x-3*y-6*z-1; 再在MATLAB命令窗口,输入命令: X1=fsolve<'fxyz',[-1,1,-1],optimset<'Display','off'>>%求第一个交点 X2=fsolve<'fxyz',[1,-1,1],optimset<'Display','off'>>%求第二个交点 例6.36求函数 在区间<-10,1>和<1,10>上的最小值点。 首先建立函数文件fx.m: functionf=f f=x-1/x+5; 上述函数文件也可用一个语句函数代替: f=inline<'x-1/x+5'> 再在MATLAB命令窗口,输入命令: fminbnd<'fx',-10,-1>%求函数在<-10,-1>内的最小值点和最小值 fminbnd 注意函数名f不用加' 例6.37设 求函数f在<0.5,0.5,0.5>附近的最小值。 建立函数文件fxyz.m: functionf=fxyz x=u<1>;y=u<2>;z=u<3>; f=x+y.^2./x/4+z.^2./y+2./z; 在MALAB命令窗口,输入命令: [U,fmin]=fminsearch<'fxyz',[0.5,0.5,0.5]>%求函数的最小值点和最小值 例6.38求解有约束最优化问题。 首先编写目标函数M文件fop.m。 functionf=fop f=0.4*x<2>+x<1>^2+x<2>^2-x<1>*x<2>+1/30*x<1>^3; 再设定约束条件,并调用fmincon函数求解此约束最优化问题。 x0=[0.5;0.5]; A=[-1,-0.5;-0.5,-1]; b=[-0.4;-0.5]; lb=[0;0]; option=optimset;option.LargeScale='off';option.Display='off'; [x,f]=fmincon<'fop',x0,A,b,[],[],lb,[],[],option> 例6.39设有初值问题: y<0>=2 试求其数值解,并与精确解相比较<精确解为y >。 <1>建立函数文件funt.m。 functiony=funt y= <2>求解微分方程。 t0=0;tf=10; y0=2; [t,y]=ode23<'funt',[t0,tf],y0>;%求数值解 y1=sqrt plot 例6.40已知一个二阶线性系统的微分方程为: 其中a=2,绘制系统的时间响应曲线和相平面图。 建立一个函数文件sys.m: functionxdot=sys xdot=[-2*x<2>;x<1>]; 取t0=0,tf=20,求微分方程的解: t0=0;tf=20; [t,x]=ode45<'sys',[t0,tf],[1,0]>; [t,x] subplot<1,2,1>;plot 2>>;%解的曲线,即t-x subplot<1,2,2>;plot 2>,x<: 1>>%相平面曲线,即x-x’ axisequal 例6.41设 将X转化为稀疏存储方式。 X=[2,0,0,0,0;0,0,0,0,0;0,0,0,5,0;0,1,0,0,-1;0,0,0,0,-5]; A=sparse 例6.42根据表示稀疏矩阵的矩阵A,产生一个稀疏存储方式矩阵B。 A=[2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12]; 例6.43求下列三对角线性代数方程组的解。 B=[1,2,0;1,4,3;2,6,1;1,6,4;0,1,2];%产生非0对角元素矩阵 d=[-1;0;1];%产生非0对角元素位置向量 A=spdiags%产生稀疏存储的系数矩阵 f=[0;3;2;1;5];%方程右边参数向量 x=%求A的全部元素的乘积。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 06 MATLAB 数值 计算 例题 源程序 汇总