MATLAB作业2参考答案(2018)

MATLAB作业二参考答案

1、试求下面线性代数方程的解析解与数值解，并检验解的正确性。

?2?10??8???5?93?2?1???1?40???3?8?4??11050??X??? ?033??2?4?63?????6?6?8?4??9?53?【求解】求出A, [A;B] 两个矩阵的秩，可见二者相同，所以方程不是矛盾方程，应该有无

穷多解。

>> A=[2,-9,3,-2,-1; 10,-1,10,5,0; 8,-2,-4,-6,3; -5,-6,-6,-8,-4]; B=[-1,-4,0; -3,-8,-4; 0,3,3; 9,-5,3]; [rank(A), rank([A B])] ans =

4 4

用下面的语句可以求出方程的解析解，并可以验证该解没有误差。 >> x0=null(sym(A));

x_analytical=sym(A)\\B; syms a; x=a*[x0 x0 x0]+x_analytical x =

[ a+967/1535, a-943/1535, a-159/1535]

[ -1535/1524*a, -1535/1524*a, -1535/1524*a]

[ -3659/1524*a-1807/1535,-3659/1524*a-257/1535,-3659/1524*a-141/1535] [ 1321/508*a+759/1535, 1321/508*a-56/1535, 1321/508*a-628/1535] [ -170/127*a-694/307, -170/127*a+719/307, -170/127*a+103/307] >> A*x-B ans =

[ 0, 0, 0] [ 0, 0, 0] [ 0, 0, 0] [ 0, 0, 0]

用数值解方法也可以求出方程的解，但会存在误差，且与任意常数a 的值有关。 >> x0=null(A); x_numerical=A\\B; syms a; x=a*[x0 x0 x0]+x_numerical; vpa(x,10) ans =

[ .2474402553*a+.1396556436, .2474402553*a-.6840666849, .2474402553*a-.1418420333]

[-.2492262414*a+.4938507789,-.2492262414*a+.7023776988e-1,-.2492262414*a+.3853511888e-1]

[ -.5940839201*a, -.5940839201*a, -.5940839201*a]

[ .6434420813*a-.7805411315, .6434420813*a-.2178190763,.6434420813*a-.5086089095]

[-.3312192394*a-1.604263460, -.3312192394*a+2.435364854, -.3312192394*a+.3867176824] >> A*x-B

ans =

[ 1/18014398509481984*a, 1/18014398509481984*a, 1/18014398509481984*a] [ -5/4503599627370496*a, -5/4503599627370496*a, -5/4503599627370496*a]

[ -25/18014398509481984*a, -25/18014398509481984*a, -25/18014398509481984*a] [ 13/18014398509481984*a, 13/18014398509481984*a, 13/18014398509481984*a]

2、求解下面的联立方程，并检验得出的高精度数值解（准解析解）的精度。

?x12?x2?1?0?22(x?2)?(x?0.5)?1?0?12

【求解】给出的方程可以由下面的语句直接求解，经检验可见，精度是相当高的。

>> [x1,x2]=solve('x1^2-x2-1=0','(x1-2)^2+(x2-0.5)^2-1=0','x1,x2') x1 =

-1.3068444845633173592407431426632-1.2136904451605911320167045558746*sqrt(-1) -1.3068444845633173592407431426632+1.2136904451605911320167045558746*sqrt(-1) 1.0673460858066897134085973128070 1.5463428833199450050728889725194 x2 =

-.76520198984055124633885565606586+3.1722093284506318231178646481143*sqrt(-1) -.76520198984055124633885565606586-3.1722093284506318231178646481143*sqrt(-1) .13922766688686144048362498805141 1.3911763127942410521940863240803

>> norm(double([x1.^2-x2-1 (x1-2).^2+(x2-0.5).^2-1])) ans =

6.058152713871457e-031

?10x1?x2?2x3?72?3、用Jacobi、Gauss-Seidel迭代法求解方程组 ??x1?10x2?2x3?83，给定初值为

??x?x?5x?423?12x(0)?(0,0,0)T。

【求解】：编写Jacobi、Gauss-Seidel函数计算，

function y=jacobi(a,b,x0)

D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1); B=D\\(L+U); f=D\\b; y=B*x0+f; n=1;

while norm(y-x0)>=1.0e-6 x0=y;

y=B*x0+f; n=n+1; end n

>> a=[10,-1,-2;-1,10,-2;-1,-1,5];b=[72;83;42]; >> jacobi(a,b,[0;0;0])

n = 17 ans =

11.0000 12.0000 13.0000

function y=seidel(a,b,x0)

D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1); G=(D-L)\\U ;f=(D-L)\\b; y=G*x0+f; n=1;

while norm(y-x0)>=1.0e-6 x0=y;

y=G*x0+f; n=n+1; end n

>> seidel(a,b,[0;0;0]) n = 10 ans =

11.0000 12.0000 13.0000

4、假设已知一组数据，试用插值方法绘制出x?(?2,4.9)区间内的光滑函数曲线，比较各

种插值算法的优劣。 -2 -1.7 -1.4 -1.1 -0.8 -0.5 -0.2 0.1 0.4 0.7 1 1.3 xi yi xi yi .10289 .11741 .13158 .14483 .15656 .16622 .17332 .1775 .17853 .17635 .17109 .16302 1.6 1.9 2.2 2.5 2.8 3.1 3.4 3.7 4 4.3 4.6 4.9 .15255 .1402 .12655 .11219 .09768 .08353 .07015 .05786 .04687 .03729 .02914 .02236 【求解】用下面的语句可以立即得出给定样本点数据的三次插值与样条插值，得出的结果如，可见，用两种插值方法对此例得出的结果几乎一致，效果均很理想。 >> 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')

5、 假设已知实测数据由下表给出，试对(x,y)在(0.1,0.1)~(1.1,1.1)区域内的点进行插值，并

用三维曲面的方式绘制出插值结果。

yi 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 x1 0.1 .83041 .83172 .83587 .84286 .85268 .86532 .88078 .89904 .92006 .94381 .97023 x2 0.2 .82727 .83249 .84345 .86013 .88251 .91049 .94396 .98276 1.0266 1.0752 1.1279 x3 0.3 .82406 .83584 .85631 .88537 .92286 .96847 1.0217 1.082 1.1482 1.2191 1.2929 x4 0.4 .82098 .84201 .87466 .91865 .97346 1.0383 1.1118 1.1922 1.2768 1.3624 1.4448 x5 0.5 .81824 .85125 .89867 .95985 1.0336 1.118 1.2102 1.3061 1.4005 1.4866 1.5564 x6 0.6 .8161 .86376 .9284 1.0086 1.1019 1.2046 1.311 1.4138 1.5034 1.5684 1.5964 x7 0.7 .81481 .87975 .96377 1.0642 1.1764 1.2937 1.4063 1.5021 1.5661 1.5821 1.5341 x8 0.8 .81463 .89935 1.0045 1.1253 1.254 1.3793 1.4859 1.5555 1.5678 1.5032 1.3473 x9 0.9 .81579 .92263 1.0502 1.1904 1.3308 1.4539 1.5377 1.5573 1.4889 1.315 1.0321 x10 1 .81853 .94959 1.1 1.257 1.4017 1.5086 1.5484 1.4915 1.3156 1.0155 .61268 x11 1.1 .82304 .9801 1.1529 1.3222 1.4605 1.5335 1.5052 1.346 1.0454 .62477 .14763 【求解】直接采用插值方法可以解决该问题，得出的插值曲面。 >> [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 命令来实现。可见，这样的插值方法更好，得出的插值 曲面更光滑。