数值分析Matlab作业.docx

1、数值分析Matlab作业数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计，电路如下： 电路中的各个电流i1,i2,i8须满足下列线性方程组：这是一个三对角方程组。设V=220V，R=27，运用追赶法，求各段电路的电流量。Matlab程序如下：function chase () %追赶法求梯形电路中各段的电流量a=input(请输入下主对角线向量a=);b=input(请输入主对角线向量b=);c=input(请输入上主对角线向量c=);d=input(请输入右端向量d=);n=input(请输入系数矩阵维数n=);u(1)=b(1);for i=2:n l(i)=a(i)/u

2、(i-1); u(i)=b(i)-c(i-1)*l(i);endy(1)=d(1);for i=2:n y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);i=n-1;while i0 x(i)=(y(i)-c(i)*x(i+1)/u(i); i=i-1;endx输入如下: chase请输入下主对角线向量a=0,-2,-2,-2,-2,-2,-2,-2;请输入主对角线向量b=2,5,5,5,5,5,5,5;请输入上主对角线向量c=-2,-2,-2,-2,-2,-2,-2,0;请输入方程组右端向量d=220/27,0,0,0,0,0,0,0;请输入系数矩阵阶数n=8运

3、行结果如下：x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477第三章14.试分别用(1)Jacobi迭代法；(2)Gauss-Seidel迭代法解线性方程组迭代初始向量。（1）雅可比迭代法程序如下：function jacobi() %Jacobi迭代法a=input(请输入系数矩阵a=);b=input(请输入右端向量b=);x0=input(请输入初始向量x0=);n=input(请输入系数矩阵阶数n=);er=input(请输入允许误差er=);N=input(请输入最大迭代次数N=);for i=1:n for j=1

4、:n if i=j d(i,j)=a(i,j); else d(i,j)=0; end endendm=eye(5)-da; %迭代矩阵g=db;x=m*x0+g;k=1;while ker x=m*x+g; k=k+1; else x return end end continueendx程序执行如下：jacobi请输入系数矩阵a=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15请输入右端向量b=12 -27 14 -17 12请输入初始向量x0=0 0 0 0 0请输入系数矩阵阶数n=5请输入允许误差er=1.0e-6

5、请输入最大容许迭代次数N=60x = 1.0000 -2.0000 3.0000 -2.00001.0000（2）高斯赛德尔迭代法程序如下：function gs_sdl() %gauss-seiddel迭代法a=input(请输入系数矩阵a=);b=input(请输入右端向量b=);x0=input(请输入初始向量x0=);n=input(请输入系数矩阵阶数n=);er=input(请输入允许误差er=);N=input(请输入最大迭代次数N=);for i=1:n for j=1:n if i=j l(i,j)=0; else l(i,j)=-a(i,j); end endendfor i

6、=1:n for j=1:n if ij u(i,j)=-a(i,j); else u(i,j)=0; end endendfor i=1:n for j=1:n if i=j d(i,j)=a(i,j); else d(i,j)=0; end endendm=(d-l)u; %迭代矩阵g=(d-l)b;x=m*x0+g;k=1;while ker x=m*x+g; k=k+1; else x return end end continueendx执行结果如下： gs_sdl请输入系数矩阵a=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3

7、 -5 -1 15请输入右端向量b=12 -27 14 -17 12请输入初始向量x0=0 0 0 0 0请输入系数矩阵阶数n=5请输入允许误差er=1.0e-6请输入最大容许迭代次数N=60x = 1.0000 -2.0000 3.0000 -2.0000 1.0000第四章已知如下矩阵，试用幂法求按模最大的特征值与特征向量。Matlab程序代码如下：function mifa ()A=input(请输入系数矩阵A=);x0=input(请输入初始列向量x0=);n=input(请输入向量维数n=);er=input(请输入允许误差er=);N=input(请输入最大容许迭代次数N=);k=

8、1;mu=0;while k mifa请输入系数矩阵A=190 66 -84 30;66 303 42 -36;336 -168 147 -112;30 -36 28 291请输入初始列向量x0=0 0 0 1请输入向量维数n=4请输入允许误差er=1.0e-6请输入最大容许迭代次数N=100lamda = 343.0000x0 = 114.3333 343.0000 -0.0000 -171.5002第五章试编写MATLAB函数实现Newton插值，要求能输出插值多项式。对函数在区间-5，5上实现10次多项式插值。Matlab程序代码如下：%此函数实现y=1/(1+4*x2)的n次Newto

9、n插值，n由调用函数时指定%函数输出为插值结果的系数向量（行向量）和插值多项式function t y=func5(n)x0=linspace(-5,5,n+1);y0=1./(1.+4.*x0.2);b=zeros(1,n+1);for i=1:n+1 s=0; for j=1:i t=1; for k=1:i if k=j t=(x0(j)-x0(k)*t; end; end; s=s+y0(j)/t; end; b(i)=s;end; t=linspace(0,0,n+1);for i=1:n s=linspace(0,0,n+1); s(n+1-i:n+1)=b(i+1).*poly(

10、x0(1:i); t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值运行结果： b Y=func5(10)b = Columns 1 through 4 -0.0000 0.0000 0.0027 -0.0000 Columns 5 through 8 -0.0514 -0.0000 0.3920 -0.0000 Columns 9 through 11 -1.1433 0.0000 1.0000Y = - (7319042784910035*x10)/147573952589676412928 + x9/18446744073709551616

11、+ (256*x8)/93425 - x7/1152921504606846976 - (28947735013693*x6)/562949953421312 - (3*x5)/72057594037927936 + (36624*x4)/93425 - (5*x3)/36028797018963968 - (5148893614132311*x2)/4503599627370496 + (7*x)/36028797018963968 + 1b为插值多项式系数向量，Y为插值多项式。插值近似值： x1=linspace(-5,5,101); x=x1(2:100); y=polyval(b,x)

12、y = Columns 1 through 12 2.7003 3.9994 4.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960 Columns 13 through 24 -0.7159 -0.7368 -0.6810 -0.5709 -0.4278 -0.2704 -0.1147 0.0270 0.1458 0.2360 0.2949 0.3227 Columns 25 through 36 0.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0

13、537 -0.0900 -0.1082 -0.1062 -0.0830 Columns 37 through 48 -0.0390 0.0245 0.1052 0.2000 0.3050 0.4158 0.5280 0.6369 0.7379 0.8269 0.9002 0.9549 Columns 49 through 60 0.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000 Columns 61 through 72 0.1052 0.0245 -0.0390 -0.0830

14、 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915 Columns 73 through 84 0.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709 Columns 85 through 96 -0.6810 -0.7368 -0.7159 -0.5960 -0.3571 0.0154 0.5274 1.1715 1.9211 2.7237 3.4926 4.0974 Columns 97 throug

15、h 994.3515 3.9994 2.7003绘制原函数和拟合多项式的图形代码：plot(x,1./(1+4.*x.2)hold allplot(x,y,r)xlabel(X)ylabel(Y)title(Runge现象)gtext(原函数)gtext(十次牛顿插值多项式)绘制结果：误差计数并绘制误差图： hold off ey=1./(1+4.*x.2)-yey = Columns 1 through 12 -2.6900 -3.9887 -4.3403 -4.0857 -3.4804 -2.7109 -1.9077 -1.1575 -0.5128 -0.0000 0.3733 0.613

16、0 Columns 13 through 24 0.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000 -0.1169 -0.2051 -0.2617 -0.2870 Columns 25 through 36 -0.2832 -0.2542 -0.2053 -0.1424 -0.0719 -0.0000 0.0674 0.1254 0.1696 0.1971 0.2062 0.1962 Columns 37 through 48 0.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0

17、.2270 -0.2379 -0.2171 -0.1649 -0.0928 Columns 49 through 60 -0.0271 0 -0.0271 -0.0928 -0.1649 -0.2171 -0.2379 -0.2270 -0.1902 -0.1349 -0.0691 0.0000 Columns 61 through 72 0.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.1254 0.0674 0.0000 -0.0719 -0.1424 Columns 73 through 84 -0.2053 -0.2542 -0.28

18、32 -0.2870 -0.2617 -0.2051 -0.1169 0.0000 0.1401 0.2943 0.4502 0.5921 Columns 85 through 96 0.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857 Columns 97 through 99 -4.3403 -3.9887 -2.6900 plot(x,ey) xlabel(X) ylabel(ey) title(Runge现象误差图)第六章16、钢包问题。炼钢唱出钢时所用的盛钢水

19、的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大经实验，钢包的容积与相应的使用次数的数据如下：使用次数x 容积y使用次数x 容积y2 106.423 108.265 109.586 109.507 109.869 110.0010 109.9311 110.5912 110.6014 110.7216 110.9017 110.7619 111.1020 111.30选用双曲线对数据进行拟合，使用最小二乘法拟合Matlab程序如下：function a=nihehanshu()x0=2 3 5 6 7 9 10 11 12 14 16 17 19 20;y0=106.4

20、2 108.26 109.58 109.50 109.86 110.00 109.93 110.59 110.60 110.72 110.90 110.76 111.10 111.30;A=zeros(2,2);B=zeros(2,1);a=zeros(2,1); x=1./x0; y=1./y0;A(1,1)=14; A(1,2)=sum(x);A(2,1)=A(1,2);A(2,2)=sum(x.2);B(1)=sum(y); B(2)=sum(x.*y);a=AB; y=1./(a(1)+a(2)*1./x0);subplot(1,2,2);plot(x0,y0-y,bd-); titl

21、e(拟合曲线误差);subplot(1,2,1);plot(x0,y0,go); hold on;x=2:0.5:20;y=1./(a(1)+a(2)*1./x);plot(x,y,r*-); legend(散点 ,拟合曲线图 1/y=a(1)+a(2)*1/x);title(最小二乘法拟合曲线);求的系数为：0.0090 0.0008则拟合曲线为 拟合曲线图、散点图、误差图如下：第七章26. 考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为 曲线关于原点对称。取a=1,参数s的变化范围-5，5，容许误差限分别是10-3和10-7。选取适当的节点个数，利用数值积分方法计算

22、曲线上点的坐标，并画出曲线的图形。程序代码如下所示：function huixuan () %用梯形公式的逐次分半算法计算回旋曲线上点的坐标er=input(请选择允许误差1.0e-3或1.0e-7：);i=1; % x向量分量的下标for s=-5:0.1:5 m=1; b=s; a=0; h=(b-a)/2; fx1=cos(a2/2); fx2=cos(b2/2); T=h*(fx1+fx2); T0=5; while abs(T-T0)3*er Fx=0; T0=T; for k=1:2(m-1) %计算新增加节点处的函数值之和 fx3=cos(a+(2*k-1)*h)2/2); Fx

23、=Fx+fx3; end T=T0/2+h*Fx; m=m+1; h=h/2; end x(i)=T; i=i+1;endj=1; %y向量分量的下标for s=-5:0.1:5 n=1; b=s; a=0; h=(b-a)/2; fy1=sin(a2/2); fy2=sin(b2/2); T=h*(fy1+fy2); T0=5; while abs(T-T0)3*er Fy=0; T0=T; for k=1:2(n-1) fy3=sin(a+(2*k-1)*h)2/2); Fy=Fy+fy3; end T=T0/2+h*Fy; n=n+1; h=h/2; end y(j)=T; j=j+1;

24、endplot(x,y,k*,x,y,k);if er=1.0e-3 title(er=1.0e-3);else title(er=1.0e-7);end程序执行结果如下： huixuan请选择允许误差1.0e-3或1.0e-7：1.0e-3 huixuan请选择允许误差1.0e-3或1.0e-7：1.0e-7 第八章20.求方程在附近的根，精确到。（1） 取，用简单迭代法计算；（2） 用加快收敛的迭代格式，计算。程序及计算过程如下：建一M文件f.m存储函数，function f=f(x)f=exp(-x);取,用简单迭代法计算，Matlab程序如下：function x,i=diedai1(

25、x0)x=f(x0);i=1;y(i)=x;while abs(x-x0)10-8 i=i+1; x0=x; x=f(x); y(i)=x; end取初始值x0=0.5,输入x,i=diedai1(0.5)得结果 x = 0.567143287611168i =30可以看出用简单收敛法经过30次迭代达到精度要求。用加速收敛法的迭代格式计算，程序如下：function x,i=diedai2(x0)w=0.625;x=w*f(x0)+(1-w)*x0;i=1;y(i)=x;while abs(x-x0)10-8 i=i+1; x0=x; x=w*f(x)+(1-w)*x; y(i)=x;end同

26、样取x0=0.5，得x = 0.567143290310401 i = 5结果比较简单迭代法和加速迭代格式的比较ix简单迭代法300.56714328761116加速的迭代格式50.567143290310401可见，加速迭代格式收敛比简单迭代格式快。第九章设有常微分方程初值问题其精确解为。选取步长使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。其中多步法需要的初值由经典RK法提供。（1）用经典四阶RK法求解，程序代码如下：function classic_rk4()n=input(请输入插值节点数n=);y(1)=1;f0(

27、1)=1; %f0=cosx+sinx为精确值h=pi/n; %步长x=0:h:pi;k=2;eps=1.0e-3;for k=1:n f0(k+1)=cos(x(k)+sin(x(k); k1=-y(k)+2*cos(x(k); k2=-(y(k)+h*k1/2)+2*cos(x(k)+h/2); k3=-(y(k)+h*k2/2)+2*cos(x(k)+h/2); k4=-(y(k)+h*k3)+2*cos(x(k)+h); y(k+1)=y(k)+h/6*(k1+2*k2+2*k3+k4); endsubplot(3,1,1);plot(x,f0,k);title(y=cosx+sinx

28、);subplot(3,1,2);plot(x,y,k);title(经典四阶RK法);subplot(3,1,3);T=y-f0; %计算经典四阶RK法的误差plot(x,T,k);title(经典四阶RK法的误差);程序执行结果如下： classic_rk4请输入插值节点数n=3000（2）用四阶Adams预测-校正算法求解，程序代码如下：function adams4()n=input(请输入插值节点数n=);h=(pi-0)/n;x=0:h:pi;for k=1:n+1 f0(k)=cos(x(k)+sin(x(k); %f0=cosx+sinx为精确值endy(1)=1;for k=2:4 %用四阶RK法获得起步值