MATLAB作业4参考答案
MATLAB课后习题集附标准答案

第2章MATLAB概论1、与其他计算机语言相比较,MA TLAB语言突出的特点是什么?答:起点高、人机界面适合科技人员、强大而简易的作图功能、智能化程度高、功能丰富,可扩展性强.2、MA TLAB系统由那些部分组成?答:开发环境、MATLAB数学函数库、MATLAB语言、图形功能、应用程序接口3、安装MATLAB时,在选择组件窗口中哪些部分必须勾选,没有勾选的部分以后如何补安装?答:在安装MATLAB时,安装内容由选择组件窗口中各复选框是否被勾选来决定,可以根据自己的需要选择安装内容,但基本平台(即MATLAB选项)必须安装.第一次安装没有选择的内容在补安装时只需按照安装的过程进行,只是在选择组件时只勾选要补装的组件或工具箱即可.4、MATLAB操作桌面有几个窗口?如何使某个窗口脱离桌面成为独立窗口?又如何将脱离出去的窗口重新放置到桌面上?答:在MATLAB操作桌面上有五个窗口,在每个窗口的右下角有两个小按钮,一个是关闭窗口的Close 按钮,一个是可以使窗口称为独立的Undock按钮,点击Undock按钮就可以使该窗口脱离桌面称为独立窗口,在独立窗口的view菜单中选择Dock,菜单项就可以将独立的窗口重新防止的桌面上.5、如何启动M文件编辑/调试器?答:在操作桌面上选择“建立新文件”或“打开文件”操作时,M文件编辑/调试器将被启动.在命令窗口中键入edit命令时也可以启动M文件编辑/调试器.6、存储在工作空间中的数组能编辑吗?如何操作?答:存储在工作空间的数组可以通过数组编辑器进行编辑:在工作空间浏览器中双击要编辑的数组名打开数组编辑器,再选中要修改的数据单元,输入修改内容即可.7、命令历史窗口除了可以观察前面键入的命令外,还有什么用途?答:命令历史窗口除了用于查询以前键入的命令外,还可以直接执行命令历史窗口中选定的内容、将选定的内容拷贝到剪贴板中、将选定内容直接拷贝到M文件中.8、如何设置当前目录和搜索路径,在当前目录上的文件和在搜索路径上的文件有什么区别?答:当前目录可以在当前目录浏览器窗口左上方的输入栏中设置,搜索路径可以通过选择操作桌面的file菜单中的Set Path菜单项来完成.在没有特别说明的情况下,只有当前目录和搜索路径上的函数和文件能够被MATLAB运行和调用,如果在当前目录上有与搜索路径上相同文件名的文件时则优先执行当前目录上的文件,如果没有特别说明,数据文件将存储在当前目录上.9、在MA TLAB中有几种获得帮助的途径?答:(1)帮助浏览器:选择view菜单中的Help菜单项或选择Help菜单中的MATLAB Help菜单项可以打开帮助浏览器.(2)help命令:在命令窗口键入“help”命令可以列出帮助主题,键入“help 函数名”可以得到指定函数的在线帮助信息.(3)lookfor命令:在命令窗口键入“lookfor 关键词”可以搜索出一系列与给定关键词相关的命令和函数.(4)模糊查询:输入命令的前几个字母,然后按Tab键,就可以列出所有以这几个字母开始的命令和函数.注意:lookfor和模糊查询查到的不是详细信息,通常还需要在确定了具体函数名称后用help命令显示详细信息.第3章 MATLAB 数值运算3.1在MA TLAB 中如何建立矩阵⎥⎦⎤⎢⎣⎡194375,并将其赋予变量a ? 答:在Command Window 窗口输入操作:>> a=[5 7 3;4 9 1]3.2有几种建立矩阵的方法?各有什么优点? 答:(1)直接输入法,如a=[1 2 3 4],优点是输入方法方便简捷;(2)通过M 文件建立矩阵,该方法适用于建立尺寸较大的矩阵,并且易于修改; (3)由函数建立,如y=sin(x),可以由MATLAB 的内部函数建立一些特殊矩阵; (4)通过数据文件建立,该方法可以调用由其他软件产生数据.3.3在进行算术运算时,数组运算和矩阵运算各有什么要求?答:进行数组运算的两个数组必须有相同的尺寸.进行矩阵运算的两个矩阵必须满足矩阵运算规则,如矩阵a 与b 相乘(a*b )时必须满足a 的列数等于b 的行数.渗釤呛俨匀谔鱉调硯錦。
MATLAB习题及参考答案经典.doc

习题:1, 计算⎥⎦⎤⎢⎣⎡=572396a 与⎥⎦⎤⎢⎣⎡=864142b 的数组乘积。
2, 对于B AX =,如果⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=753467294A ,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=282637B ,求解X 。
3, 已知:⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321a ,分别计算a 的数组平方和矩阵平方,并观察其结果。
4, 角度[]604530=x ,求x 的正弦、余弦、正切和余切。
(应用sin,cos,tan.cot)5, 将矩阵⎥⎦⎤⎢⎣⎡=7524a 、⎥⎦⎤⎢⎣⎡=3817b 和⎥⎦⎤⎢⎣⎡=2695c 组合成两个新矩阵: (1)组合成一个4⨯3的矩阵,第一列为按列顺序排列的a 矩阵元素,第二列为按列顺序排列的b 矩阵元素,第三列为按列顺序排列的c 矩阵元素,即 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡237912685574(2)按照a 、b 、c 的列顺序组合成一个行矢量,即 []2965318772546, 将(x -6)(x -3)(x -8)展开为系数多项式的形式。
(应用poly,polyvalm)7, 求解多项式x 3-7x 2+2x +40的根。
(应用roots)8, 求解在x =8时多项式(x -1)(x -2) (x -3)(x -4)的值。
(应用poly,polyvalm)9, 计算多项式9514124234++--x x x x 的微分和积分。
(应用polyder,polyint ,poly2sym)10, 解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡66136221143092x 。
(应用x=a\b)11, 求欠定方程组⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡5865394742x 的最小范数解。
(应用pinv) 12, 矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=943457624a ,计算a 的行列式和逆矩阵。
(应用det,inv)13, y =sin(x ),x 从0到2π,∆x =0.02π,求y 的最大值、最小值、均值和标准差。
matlab课后习题答案第四章

第4章数值运算习题 4 及解答1 根据题给的模拟实际测量数据的一组t和)(t y试用数值差分diff或数值梯度gradient指令计算)(t y',然后把)(t y和)(t y'曲线绘制在同一张图上,观察数值求导的后果。
(模拟数据从prob_data401.mat获得)〖目的〗●强调:要非常慎用数值导数计算。
●练习mat数据文件中数据的获取。
●实验数据求导的后果●把两条曲线绘制在同一图上的一种方法。
〖解答〗(1)从数据文件获得数据的指令假如prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat(2)用diff求导的指令dt=t(2)-t(1);yc=diff(y)/dt; %注意yc的长度将比y短1plot(t,y,'b',t(2:end),yc,'r')(3)用gradent 求导的指令(图形与上相似)dt=t(2)-t(1);yc=gradient(y)/dt;plot(t,y,'b',t,yc,'r')grid on〖说明〗● 不到万不得已,不要进行数值求导。
● 假若一定要计算数值导数,自变量增量dt 要取得比原有数据相对误差高1、2个量级以上。
● 求导会使数据中原有的噪声放大。
2 采用数值计算方法,画出dt tt x y x ⎰=0sin )(在]10 ,0[区间曲线,并计算)5.4(y 。
〖提示〗● 指定区间内的积分函数可用cumtrapz 指令给出。
● )5.4(y 在计算要求不太高的地方可用find 指令算得。
〖目的〗● 指定区间内的积分函数的数值计算法和cumtrapz 指令。
● find 指令的应用。
〖解答〗dt=1e-4;t=0:dt:10;t=t+(t==0)*eps;f=sin(t)./t;s=cumtrapz(f)*dt;plot(t,s,'LineWidth',3)ii=find(t==4.5);s45=s(ii)s45 =1.65413 求函数x ex f 3sin )(=的数值积分⎰=π0 )(dx x f s ,并请采用符号计算尝试复算。
matlab习题四答案

习题四1、从键盘输入一个4位数,按如下规则加密输出,加密规则:每位数字都加7,然后用和除以10的余数来代替概数,一三位交换,二四位交换。
function[a1 a2 a3 a4]=jiami(s)s=input('输入需要加密的数:');x1=rem(s,10);x2=fix(rem(s,100)/10);x3=fix(rem(s,1000)/100);x4=fix(s/1000);x1x2x3x4a1=rem(x3+7,10);a2=rem(x4+7,10);a3=rem(x1+7,10);a4=rem(x2+7,10);a1a2a3a4运行:>> jiami输入需要加密的数:4567x1 =7x2 =6x3 =5x4 =4a1 =2a2 =1a3 =4a4 =32、分别用if和switch语句实现以下计算,其中a,b,c的值是从键盘输入。
(1)用if函数;function y=fangcheng(a,b,c,x)a=input('inputa');b=input('inputb');c=input('inputc');x=input('enterX:');if (x>=0.5&x<1.5)y=a*x*x+b*x+c;elseif (x>=1.5&x<3.5)y=a*sin(x)+x;elseif (x>=3.5&x<5.5)y=log10(abs(b+c/x))elsedisp('error');end> fangchenginputa5inputb6inputc7enterX:2.4ans =5.7773(2)用switch函数;function y=fangcheng2(a,b,c,x)a=input('inputa');b=input('inputb');c=input('inputc');x=input('enterX:');switch fix(x/2)case 0y=a*x*x+b*x+c;case 1y=a*sin(x)+x;case 2y=log10(abs(b+c/x));otherwise>> fangcheng2inputa5inputb8inputc66enterX:5.3ans =1.31083、产生20个随机正数,输出其中小于平均值的偶数;a=randint(1,20,[30,100]);sum=0;for i=1:20sum=sum+a(i);endj=sum/20;for i=1:20;if (rem(a(i),2)==0&a(i)<j)a(i)endend>> qiuheans =42ans =30ans =44ans =50ans =344、输入20个随机整数,求出其中最大和最小的数。
电机与拖动基础及MATLAB仿真习题答案(第四章)

4-14 一台直流电动机技术数据如下:额定功率PN=40kW ,额定电压UN=220V ,额定转速nN=1500r/min ,额定效率η=%,求电动机的额定电流和额定负载时的输入功率 解:(1)额定电流(2)输入功率4-15 一台直流发电机技术数据如下:额定功率PN=82kW ,额定电压UN=230V ,额定转速nN=970r/min ,额定效率η=90%,求发电机的额定电流和额定负载时的输入功率 解:(1)额定电流(2)输入功率4-16 已知一台直流电机极对数p=2,槽数Z 和换向片数K 均等于22,采用单叠绕组。
试求:(1)绕组各节距;(2)并联支路数。
解:(1)第一节距5424222y 1=-=±=εp z ,为短距绕组。
单叠绕组的合成节距及换向器节距均为1,即1y ==k y第二节距415y 12=-=-=y y(2)并联支路数等于磁极数,为4。
4-17 已知直流电机极数2p=6,电枢绕组总导体数N=400,电枢电流Ia=10A ,气隙每极磁通Φ=×10-2Wb ,试求:(1)采用单叠绕组时电枢所受电磁转矩;(2)绕组改为单波保持支路电流ia 不变时的电磁转矩。
解: 电枢绕组为单叠绕组时,并联支路对数a=p=3,电磁转矩 m N I a pN T a ⋅=⨯⨯⨯⨯⨯=Φ=38.1310021.0314.3240032π 如果把电枢绕组改为单波绕组, 保持支路电流a i 的数值不变,则电磁转矩也不变,仍A U P I N N N N 79.207875.022010403=⨯⨯==ηkWI U P N N 71.4579.2072201=⨯=⨯=A U P I N N N 5.35623010823=⨯==KW P P N 11.911==η为m N ⋅,因为无论是叠绕组还是波绕组,所有导体产生的电磁转矩的方向是一致的,保持支路电流a i 不变,就保持了导体电流不变,也就保持了电磁转矩不变。
matlab第四章课后作业解答

matlab第四章课后作业解答第四章习题解答1、求下列多项式的所有根,并进行验算。
(3)267235865x x x x-+-(4)4)32(3-+x 解:>> p=zeros(1,24); >> p(1)=5;p(17)=-6;p(18)=8;p(22)=-5; >> root=roots(p)root =0.97680.9388 + 0.2682i0.9388 - 0.2682i0.8554 + 0.5363i0.8554 - 0.5363i0.6615 + 0.8064i0.6615 - 0.8064i0.3516 + 0.9878i0.3516 - 0.9878i-0.0345 + 1.0150i-0.0345 - 1.0150i-0.4609 + 0.9458i-0.4609 - 0.9458i-0.1150 + 0.8340i-0.1150 - 0.8340i-0.7821 + 0.7376i-0.7821 - 0.7376i-0.9859 + 0.4106i-0.9859 - 0.4106i-1.0416-0.7927>> polyval(p,root)ans =1.0e-012 *-0.07120.0459 - 0.0081i0.0459 + 0.0081i-0.0419 + 0.0444i-0.0419 - 0.0444i0.0509 + 0.0929i0.0509 - 0.0929i-0.2059 + 0.0009i-0.2059 - 0.0009i-0.0340 + 0.0145i-0.0340 - 0.0145i0.1342 + 0.0910i0.1342 - 0.0910i0.0025 + 0.0027i0.0025 - 0.0027i-0.0077 + 0.4643i-0.0077 - 0.4643i-0.3548 - 0.1466i-0.3548 + 0.1466i-0.0251-0.0073(4) >> p1=[2 3];>> p=conv(conv(p1,p1),p1)-[0 0 0 4]; >> root=roots(p)root =-1.8969 + 0.6874i-1.8969 - 0.6874i-0.7063>> polyval(p,root)ans =1.0e-014 *-0.7105 - 0.6217i-0.7105 + 0.6217i6、求解下列方程组在区域1,0<<βα内的解-=+=.sin 2.0cos 7.0,cos 2.0sin 7.0βαββαα 解:以初值)5.0,5.0(),(00=βα进行求解>> fun=inline('[0.7*sin(x(1))+0.2*cos(x(2))-x(1),0.7*cos(x(1))-0.2*sin(x(2))-x(2)]');>> [x,f,h]=fsolve(fun,[0.5 0.5])Optimization terminated: first-order optimality is less than options.TolFun.x =0.5265 0.5079f =1.0e-007 *-0.1680 -0.2712h =1因而,该方程组的近似根为5079.0,5265.0==βα。
matlab课后习题与答案

习题二1.如何理解“矩阵是MATLAB最基本的数据对象”?答:因为向量可以看成是仅有一行或一列的矩阵,单个数据(标量)可以看成是仅含一个元素的矩阵,故向量和单个数据都可以作为矩阵的特例来处理。
因此,矩阵是MATLAB最基本、最重要的数据对象。
2.设A和B是两个同维同大小的矩阵,问:(1)A*B和A.*B的值是否相等?答:不相等。
(2)A./B和B.\A的值是否相等?答:相等。
(3)A/B和B\A的值是否相等?答:不相等。
(4)A/B和B\A所代表的数学含义是什么?答:A/B等效于B的逆右乘A矩阵,即A*inv(B),而B\A等效于B矩阵的逆左乘A矩阵,即inv(B)*A。
3.写出完成下列操作的命令。
(1)将矩阵A第2~5行中第1, 3, 5列元素赋给矩阵B。
答:B=A(2:5,1:2:5); 或B=A(2:5,[1 3 5])(2)删除矩阵A的第7号元素。
答:A(7)=[](3)将矩阵A的每个元素值加30。
答:A=A+30;(4)求矩阵A的大小和维数。
答:size(A);ndims(A);(5)将向量t的0元素用机器零来代替。
答:t(find(t==0))=eps;(6)将含有12个元素的向量x转换成34矩阵。
答:reshape(x,3,4);(7)求一个字符串的ASCII码。
答:abs(‘123’); 或double(‘123’);(8) 求一个ASCII 码所对应的字符。
答:char(49);4. 下列命令执行后,L1、L2、L3、L4的值分别是多少?A=1:9;B=10-A;...L1=A==B;L2=A<=5;L3=A>3&A<7;L4=find(A>3&A<7);答:L1的值为[0, 0, 0, 0, 1, 0, 0, 0, 0]L2的值为[1, 1, 1, 1, 1, 0, 0, 0, 0]L3的值为[0, 0, 0, 1, 1, 1, 0, 0, 0]L4的值为[4, 5, 6]5. 已知23100.7780414565532503269.5454 3.14A -⎡⎤⎢⎥-⎢⎥=⎢⎥⎢⎥-⎣⎦完成下列操作:(1) 取出A 的前3行构成矩阵B ,前两列构成矩阵C ,右下角32⨯子矩阵构成矩阵D ,B 与C 的乘积构成矩阵E 。
matlab实验4_循环结构程序设计_参考解答

高教社刘卫国《MATLAB 程序设计与应用》(第二版)实验参考解答
实验四:循环结构程序设计
clear all; close all; clc; y=0; n=1; while y<3
y=y+1/(2*n-1); n=n+1; end n=n-1; y=y-1/(2*n-1) n=n-1
运行结果:
y= 2.9944
x1=a/(b+x0); y=abs(x1-x0); x0=x1;
MATLAB 语言
3
WangChengyou © 山东大学威海分校
高教社刘卫国《MATLAB 程序设计与应用》(第二版)实验参考解答
实验四:循环结构程序设计
n=n+1; end x1 r1=(-b+sqrt(b^2+4*a))/2 r2=(-b-sqrt(b^2+4*a))/2
clear all; close all; clc; prime=[]; for n=2:49
x=n*(n+1)-1; signx=0; for k=2:sqrt(x)
MATLAB 语言
5
WangChengyou © 山东大学威海分校
高教社刘卫国《MATLAB 程序设计与应用》(第二版)实验参考解答
根据
π2 6
=1 12
+
1 22
+1 32
+"+
1 n2
,求 π 的近似值。当
n
分别取
100、1000、10000
时,
结果是多少?
要求:分别用循环结构和向量运算(使用 sum 函数)来实现。 答:程序设计:
(1) 循环结构
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
MATLAB 作业四参考答案1、 用2sin(103)y t =+在(0,3)区间内生成一组较稀疏的数据,并用一维数据插值的方法对给出的数据进行曲线拟合,并将结果与理论曲线相比较。
【求解】类似于上面的例子,可以用几乎一致的语句得出样本数据和插值效果。
>> t=0:0.2:3;y=sin(10*t.^2+3); plot(t,y,'o')ezplot('sin(10*t^2+3)',[0,3]); hold on x1=0:0.001:3; y1=interp1(t,y,x1,'spline'); plot(x1,y1)由于曲线本身变换太大,所以在目前选定的样本点下是不可能得出理想插值效果的,因为样 本数据提供的信息量不够。
为了得到好的插值效果,必须增大样本数据的信息量,对本例来 说,必须在快变化区域减小样本点的步长。
>> hold offt=[0:0.1:1,1.1:0.04:3]; y=sin(10*t.^2+3); plot(t,y,'o') ezplot('sin(10*t^2+3)',[0,3]); hold on x1=0:0.001:3; y1=interp1(t,y,x1,'spline'); plot(x1,y1)2、 用242231(,)sin()3x y f x y e xy x y x y--=++原型函数生成一组网络数据或随机数据,分别拟合出曲面,并和原曲面进行比较。
【求解】由下面的语句可以直接生成一组网格数据,用下面语句还可以还绘制出给定样本点是三维表面图。
>> [x,y]=meshgrid(0.2:0.2:2);z=exp(-x.^2-y.^4).*sin(x.*y.^2+x.^2.*y)./(3*x.^3+y); surf(x,y,z)选择新的密集网格,则可以通过二元插值得出插值曲面。
对比插值结果和新网格下的函数值精确解,则可以绘制出绝对插值误差曲面。
由插值结果可见精度是令人满意的。
>> [x1,y1]=meshgrid(0.2:0.02:2); z1=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1)>> z0=exp(-x1.^2-y1.^4).*sin(x1.*y1.^2+x1.^2.*y1)./(3*x1.^3+y1); surf(x1,y1,abs(z1-z0))现在假设已知的样本点不是网格形式分布的,而是随机分布的,则可以用下面语句生成样本点,得出分布的二维、三维示意图。
>> x=0.2+1.8*rand(400,1); y=0.2+1.8*rand(400,1); % 仍生成(0.2,2) 区间的均匀分布随机数z=exp(-x.^2-y.^4).*sin(x.*y.^2+x.^2.*y)./(3*x.^3+y); plot(x,y,'x')figure, plot3(x,y,z,'x')利用下面的语句可以得出三维插值结果,同时可以绘制出插值的绝对误差曲面,可见插值结果还是很好的,但由于边界样本点信息不能保证,所以不能像网格数据那样对(0.2,2) 区域,而只能选择(0.3,1.9) 区域进行插值。
>> [x1,y1]=meshgrid(0.3:0.02:1.9);z1=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z1)>> z0=exp(-x1.^2-y1.^4).*sin(x1.*y1.^2+x1.^2.*y1)./(3*x1.^3+y1);surf(x1,y1,abs(z1-z0))x∈-区间内的光滑函数曲线,比较各种插值算法3、假设已知一组数据,试用插值方法绘制出(2,4.9)【求解】用下面的语句可以立即得出给定样本点数据的三次插值与样条插值,得出的结果如,可见,用两种插值方法对此例得出的结果几乎一致,效果均很理想。
>> x=[-2,-1.7,-1.4,-1.1,-0.8,-0.5,-0.2,0.1,0.4,0.7,1,1.3,...1.6,1.9,2.2,2.5,2.8,3.1,3.4,3.7,4,4.3,4.6,4.9];y=[0.10289,0.11741,0.13158,0.14483,0.15656,0.16622,0.17332,...0.1775,0.17853,0.17635,0.17109,0.16302,0.15255,0.1402,...0.12655,0.11219,0.09768,0.08353,0.07019,0.05786,0.04687,...0.03729,0.02914,0.02236];x0=-2:0.02:4.9;y1=interp1(x,y,x0,'cubic');y2=interp1(x,y,x0,'spline');plot(x0,y1,':',x0,y2,x,y,'o')x y在(0.1,0.1)~(1.1,1.1)区域内的点进行插值,并用三维曲面4、假设已知实测数据由下表给出,试对(,)>> [x,y]=meshgrid(0.1:0.1:1.1);z=[0.83041,0.82727,0.82406,0.82098,0.81824,0.8161,0.81481,0.81463,0.81579,0.81853,0.82304 ;0.83172,0.83249,0.83584,0.84201,0.85125,0.86376,0.87975,0.89935,0.92263,0.94959,0.9801; 0.83587,0.84345,0.85631,0.87466,0.89867,0.9284,0.96377,1.0045,1.0502,1.1,1.1529;0.84286,0.86013,0.88537,0.91865,0.95985,1.0086,1.0642,1.1253,1.1904,1.257,1.3222;0.85268,0.88251,0.92286,0.97346,1.0336,1.1019,1.1764,1.254,1.3308,1.4017,1.4605;0.86532,0.91049,0.96847,1.0383,1.118,1.2046,1.2937,1.3793,1.4539,1.5086,1.5335;0.88078,0.94396,1.0217,1.1118,1.2102,1.311,1.4063,1.4859,1.5377,1.5484,1.5052;0.89904,0.98276,1.082,1.1922,1.3061,1.4138,1.5021,1.5555,1.5573,1.4915,1.346;0.92006,1.0266,1.1482,1.2768,1.4005,1.5034,1.5661,1.5678,1.4889,1.3156,1.0454;0.94381,1.0752,1.2191,1.3624,1.4866,1.5684,1.5821,1.5032,1.315,1.0155,0.62477;0.97023,1.1279,1.2929,1.4448,1.5564,1.5964,1.5341,1.3473,1.0321,0.61268,0.14763];[x1,y1]=meshgrid(0.1:0.02:1.1);z1=interp2(x,y,z,x1,y1,'spline');surf(x1,y1,z1)axis([0.1,1.1,0.1,1.1,min(z1(:)),max(z1(:))])其实,若光需要插值曲面而不追求插值数值的话,完全可以直接采用MATLAB 下的shading interp 命令来实现。
可见,这样的插值方法更好,得出的插值曲面更光滑。
>> surf(x,y,z); shading interp5、习题3和4给出的数据分别为一元数据和二元数据,试用分段三次样条函数和B样条函数对其进行拟合。
【求解】先考虑习题4,相应的三次样条插值和B-样条插值原函数与导数函数分别为:>> x=[-2,-1.7,-1.4,-1.1,-0.8,-0.5,-0.2,0.1,0.4,0.7,1,1.3,...1.6,1.9,2.2,2.5,2.8,3.1,3.4,3.7,4,4.3,4.6,4.9];y=[0.10289,0.11741,0.13158,0.14483,0.15656,0.16622,0.17332,...0.1775,0.17853,0.17635,0.17109,0.16302,0.15255,0.1402,...0.12655,0.11219,0.09768,0.08353,0.07019,0.05786,0.04687,...0.03729,0.02914,0.02236];S=csapi(x,y); S1=spapi(6,x,y);fnplt(S); hold on; fnplt(S1)>> Sd=fnder(S); Sd1=fnder(S1);fnplt(Sd), hold on; fnplt(Sd1)再考虑习题5 中的数据,原始数据不能直接用于样条处理,因为meshgrid() 函数产生的数据格式与要求的ndgrid() 函数不一致,所以需要对数据进行处理,其中需要的x 和y 均应该是向量,而z 是原来z矩阵的转置,所以用下面的语句可以建立起三次样条和B-样条的插值模型,函数的表面图所示,可见二者得出的结果很接近。
>> [x,y]=meshgrid(0:0.1:1.1);z=[0.83041,0.82727,0.82406,0.82098,0.81824,0.8161,0.81481,0.81463,0.81579,0.81853,0.82304;0.83172,0.83249,0.83584,0.84201,0.85125,0.86376,0.87975,0.89935,0.92263,0.94959,0.9801; 0.83587,0.84345,0.85631,0.87466,0.89867,0.9284,0.96377,1.0045,1.0502,1.1,1.1529; 0.84286,0.86013,0.88537,0.91865,0.95985,1.0086,1.0642,1.1253,1.1904,1.257,1.3222; 0.85268,0.88251,0.92286,0.97346,1.0336,1.1019,1.1764,1.254,1.3308,1.4017,1.4605; 0.86532,0.91049,0.96847,1.0383,1.118,1.2046,1.2937,1.3793,1.4539,1.5086,1.5335; 0.88078,0.94396,1.0217,1.1118,1.2102,1.311,1.4063,1.4859,1.5377,1.5484,1.5052; 0.89904,0.98276,1.082,1.1922,1.3061,1.4138,1.5021,1.5555,1.5573,1.4915,1.346; 0.92006,1.0266,1.1482,1.2768,1.4005,1.5034,1.5661,1.5678,1.4889,1.3156,1.0454; 0.94381,1.0752,1.2191,1.3624,1.4866,1.5684,1.5821,1.5032,1.315,1.0155,0.62477; 0.97023,1.1279,1.2929,1.4448,1.5564,1.5964,1.5341,1.3473,1.0321,0.61268,0.14763]; >> x0=[0.0:0.1:1]; y0=x0; z=z'; S=csapi({x0,y0},z); fnplt(S)figure; S1=spapi({5,5},{x0,y0},z); fnplt(S1) >> S1x=fnder(S1,[0,1]); fnplt(S1x)figure; S1y=fnder(S1,[0,1]); fnplt(S1y)6、 重新考虑习题3中给出的数据,试考虑用多项式插值的方法对其数据进行逼近,并选择一个能较好拟合原数据的多项式阶次。