数值计算方法实践作业.docx

1、 用Gauss-Seidel迭代法求解方程组。1.4449x1+0.7948x2+0.8801x3=10.6946x1+1.9568x2+0.1730x3= 00.6213x1+0.5526x2+1.9797x3=1解: function x,n=gauseidel(A,b,x0,eps,M)if nargin=3eps= 1.0e-6;M = 200;elseif nargin = 4M = 200;elseif nargin=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp(Warning: );return;endend A=1.446 0.7948 0.8801;0.6946 1.9568 0.1730;0.6213 0.5526 1.9797; b=1 0 1; x0=zeros(3,1); x,n=gauseidel(A,b,x0)x = 0.5891-0.24340.3882n =112、 用Jacobi迭代法求解方程组。0.9889x1+0.0005x2-0.0002x3=1-0.0046x1+0.9946x2+0.0077= 0-0.0002x1+0.0092x2+0.9941x3=1解: function x,n=jacobi(A,b,x0,eps,varargin)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(Warning: ); return; endend A=0.9989 -0.0005 -0.0002;-0.0046 0.9946 0.0077;-0.0002 0.0092 0.9941; b=1 0 1; x0=zeros(3,1); x,n=jacobil(A,b,x0)x = 1.0114-0.00311.0062n =43、 求出下表所列各数据的拉格朗日插值多项式，并计算x=1.6时y的值。X11.21.82.54y0.84150.93200.97380.5985-0.7568解：function f = Language(x,y,x0)syms t;if(length(x) = length(y) n = length(x); else disp(xy); return;end % f = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); % end; f = f + l; % simplify(f); % if(i=n) if(nargin = 3) f = subs(f,t,x0); % else f = collect(f); % f = vpa(f,6); %6 end endend x=1 1.2 1.8 2.5 4; y=0.8415 0.9320 0.9738 0.5985 -0.7568; f=Language(x,y) f = 0.0328112*t4 - 0.204917*t3 - 0.0145485*t2 + 1.05427*t - 0.0261189 f=Language(x,y,1.6)f = 0.99924、 求出下表所列各数据的均差形式的牛顿插值多项式，并计算x=1.6时y的值。X11.21.82.54y11.443.246.2516解：function f = Newton(x,y,x0)syms t; if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(xy); return;end f = y(1);y1 = 0;l = 1; for(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f); y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,x0); else f = collect(f); % f = vpa(f, 6); end endend x=1 1.2 1.8 2.5 4; y=1 1.44 3.24 6.25 16; f=Newton(x,y) f = 2.11471e-16*t4 - 1.69177e-15*t3 + 1.0*t2 - 4.82154e-15*t + 1.82711e-15 f=Newton(x,y,1.6)f = 2.56005、 求出下表所列各数据的前向差分形式的牛顿插值多项式，并计算x=1.55时y的值。X11.21.41.61.8y0.84150.93200.98540.99960.9738解：function f = Newtonforward(x,y,x0)syms t; if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(xy); return;end f = y(1);y1 = 0; xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp(); return;end for(i=1:n-1) for(j=1:n-i) y1(j) = y(j+1)-y(j); end c(i) = y1(1); l = t; for(k=1:i-1) l = l*(t-k); end; f = f + c(i)*l/factorial(i); simplify(f); y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,(x0-x(1)/(x(2)-x(1); else f = collect(f); f = vpa(f, 6); end endend x=1:0.2:1.8; y=0.8415 0.9320 0.9854 0.9996 0.9738; f=Newtonforward(x,y) f=Newtonforward(x,y) f = 0.0000541667*t4 - 0.000675*t3 - 0.0169042*t2 + 0.108025*t + 0.8415 f=Newtonforward(x,y,1.55)f = 0.99986、 求出下表所列各数据的高斯插值多项式，并计算x=1.5时y的值。X11.21.41.61.8y0.84150.93200.98540.99960.9738解：function f = Gauss(x,y,x0) if(length(x) = length(y) n = length(x);else disp(xy); return;end xx =linspace(x(1),x(n),(x(2)-x(1); if( mod(n,2) =1) if(nargin = 2) f = GStirling(x,y,n); else if(nargin = 3) f = GStirling(x,y,n,x0); end endelse if(nargin = 2) f = GBessel(x,y,n); else if(nargin = 3) f = GBessel(x,y,n,x0); end endend function f = GStirling(x,y,n,x0)syms t;nn = (n+1)/2;f = y(nn); for(i=1:n-1) for(j=i+1:n) y1(j) = y(j)-y(j-1); end if(mod(i,2)=1) c(i) = (y1(i+n)/2)+y1(i+n+2)/2)/2; else c(i) = y1(i+n+1)/2)/2; end if(mod(i,2)=1) l = t+(i-1)/2; for(k=1:i-1) l = l*(t+(i-1)/2-k); end else l_1 = t+i/2-1; l_2 = t+i/2; for(k=1:i-1) l_1 = l_1*(t+i/2-1-k); l_2 = l_2*(t+i/2-k); end l = l_1 + l_2; end l = l/factorial(i); f = f + c(i)*l; simplify(f); f = vpa(f, 6); y = y1; if(i=n-1) if(nargin = 4) f = subs(f,t,(x0-x(nn)/(x(2)-x(1); end endend function f = GBessel(x,y,n,x0)syms t;nn = n/2;f = (y(nn)+y(nn+1)/2; for(i=1:n-1) for(j=i+1:n) y1(j) = y(j)-y(j-1); end if(mod(i,2)=1) c(i) = y1(i+n+1)/2)/2; else c(i) = (y1(i+n)/2)+y1(i+n+2)/2)/2; end if(mod(i,2)=0) l = t+i/2-1; for(k=1:i-1) l = l*(t+i/2-1-k); end else l_1 = t+(i-1)/2; l_2 = t+(i-1)/2-1; for(k=1:i-1) l_1 = l_1*(t+(i-1)/2-k); l_2 = l_2*(t+(i-1)/2-1-k); end l = l_1 + l_2; end l = l/factorial(i); f = f + c(i)*l; simplify(f); f = vpa(f, 6); y = y1; if(i=n-1) if(nargin = 4) f = subs(f,t,(x0-x(nn)/(x(2)-x(1); end endend x=1:0.2:1.8; y=0 0.2630 0.4854 0.6781 0.8480; f=Gauss(x,y) f=Gauss(x,y) f = 0.20755*t - 0.