当前位置:文档之家› 数值分析Matlab作业

数值分析Matlab作业

数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计,电路如下:电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组:121232343454565676787822/252025202520252025202520250i i V Ri i ii i ii i ii i ii i ii i ii i-=-+-=-+-=-+-=-+-=-+-=-+-=-+=这是一个三对角方程组。

设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:nl(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*l(i);endy(1)=d(1);for i=2:ny(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);i=n-1;while i>0x(i)=(y(i)-c(i)*x(i+1))/u(i);i=i-1;endx输入如下:请输入下主对角线向量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 运行结果如下:x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477第三章14.试分别用(1)Jacobi 迭代法;(2)Gauss-Seidel 迭代法解线性方程组1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦ 迭代初始向量(0)(0,0,0,0,0)T x =。

(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:n if i==jd(i,j)=a(i,j); elsed(i,j)=0; end end endm=eye(5)-d\a; %迭代矩阵 g=d\b;x=m*x0+g; k=1;while k<=N %进行迭代 for i=1:5if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;xreturnendendcontinueendx程序执行如下:>>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请输入最大容许迭代次数N=60x =1.0000-2.00003.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:nfor j=1:nif i<=jl(i,j)=0;elsel(i,j)=-a(i,j);endendendfor i=1:nfor j=1:nif i<ju(i,j)=-a(i,j);elseu(i,j)=0;endendendfor i=1:nfor j=1:nif i==jd(i,j)=a(i,j);elsed(i,j)=0;endendendm=(d-l)\u; %迭代矩阵g=(d-l)\b;x=m*x0+g;k=1;while k<=Nfor i=1:5if max(abs(x(i)-x0(i))) >erx=m*x+g;k=k+1;elsexreturnendendcontinueendx执行结果如下:>> 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 -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.00003.0000-2.00001.0000已知如下矩阵,试用幂法求按模最大的特征值与特征向量。

190668430663034236336168147112303628291-⎡⎤⎢⎥-⎢⎥⎢⎥--⎢⎥-⎣⎦Matlab 程序代码如下:function mifa ()A=input('请输入系数矩阵A='); x0=input('请输入初始列向量x0='); n=input('请输入向量维数n='); er=input('请输入允许误差er=');N=input('请输入最大容许迭代次数N='); k=1; mu=0;while k<=N for t=1:nif abs(x0(t))==max(abs(x0)) alfa=x0(t);xb=t; %最大的x0(i )的下标 end endy=x0./alfa; x0=A*y;lamda=x0(xb); k=k+1; endlamda %按模最大的特征值x0 %按模最大的特征值对应的特征向量 程序执行结果如下: >> 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=100 lamda = 343.0000 x0 =114.3333 343.0000 -0.0000 -171.5002试编写MATLAB函数实现Newton插值,要求能输出插值多项式。

对函数在区间[-5,5]上实现10次多项式插值。

Matlab程序代码如下:%此函数实现y=1/(1+4*x^2)的n次Newton插值,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+1s=0;for j=1:it=1;for k=1:iif k~=jt=(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:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(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.0000Columns 5 through 8-0.0514 -0.0000 0.3920 -0.0000Columns 9 through 11-1.1433 0.0000 1.0000Y =- (7319042784910035*x^10)/147573952589676412928 + x^9/18446744073709551616 + (256*x^8)/93425 - x^7/1152921504606846976 - (28947735013693*x^6)/562949953421312 - (3*x^5)/72057594037927936 + (36624*x^4)/93425 - (5*x^3)/36028797018963968 - (5148893614132311*x^2)/4503599627370496 + (7*x)/36028797018963968 + 1b为插值多项式系数向量,Y为插值多项式。

插值近似值:x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y =Columns 1 through 122.70033.99944.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960Columns 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.3227Columns 25 through 360.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830Columns 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.9549Columns 49 through 600.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000Columns 61 through 720.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915Columns 73 through 840.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709Columns 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.0974Columns 97 through 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 offey=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.6130Columns 13 through 240.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000-0.1169 -0.2051 -0.2617 -0.2870Columns 25 through 36-0.2832 -0.2542 -0.2053 -0.1424 -0.0719 -0.0000 0.0674 0.12540.1696 0.1971 0.2062 0.1962Columns 37 through 480.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270-0.2379 -0.2171 -0.1649 -0.0928Columns 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.0000Columns 61 through 720.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.12540.0674 0.0000 -0.0719 -0.1424Columns 73 through 84-0.2053 -0.2542 -0.2832 -0.2870 -0.2617 -0.2051 -0.1169 0.00000.1401 0.2943 0.4502 0.5921 Columns 85 through 960.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.6900plot(x,ey) xlabel('X') ylabel('ey')title('Runge 现象误差图')第六章16、钢包问题。

相关主题