数值分析作业-matlab上机作业

1、数值分析Matlab上机作业学院：班级：老师：姓名：学号：第二章 解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a，b，c，d输出结果x，均为数组形式function x=Zhuiganfa(a,b,c,d)%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。%定义三对角矩阵A的各组成单元。方程为Ax=d%b为A的对角线元素(1n)，a为-1对角线元素(2n)，c为+1对角线元素(1n-1)。% A=2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5% a=-1 -2 -3;c=-1 -2 -3;b=2 3 4 5

2、;d=6 1 -2 1;n=length(b);u(1)=b(1);y(1)=d(1);for i=2:n l(i)=a(i-1)/u(i-1);%先求l(i) u(i)=b(i)-c(i-1)*l(i);%再求u(i) %A=LU,Ax=LUx=d,y=Ux, %Ly=d,由于L是下三角矩阵，对角线均为1，所以可求y(i) y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=(n-1):-1:1 %Ux=y,由于U是上三角矩阵，所以可求x(i) x(i)=(y(i)-c(i)*x(i+1)/u(i);end2、输入已知参数a=2 2 2 2 2 2 2

3、;b=2 5 5 5 5 5 5 5;c=2 2 2 2 2 2 2;d=220/27 0 0 0 0 0 0 0;3、按定义格式调用函数x=zhuiganfa(a,b,c,d)4、输出结果x=8.9105 -4.5030 2.8471 -1.1148 0.9400 -0.7350 0.3976 -0.7591第15题【解】1、编写一个程序生成题目条件生成线性方程组Ax=b的系数矩阵A和右端项量b，分别定义矩阵A、B、a、b分别表示系数矩阵，其中或分别构成A、B对应右端项量分别a、b。程序如下：clear,clc;n=5;%定义A矩阵A=zeros(n,n);for i=1:n x=1+0.1

4、*i; for j=1:n A(i,j)=x(j-1); endend%定义B矩阵B=zeros(n,n);for i=1:n for j=1:n B(i,j)=1/(i+j-1); endend%定义a向量，其中Ax=afor i=1:n x=1+0.1*i; a(i)=0; for j=1:n a(i)=x(j-1)+a(i); endend %定义b向量，其中Bx=bfor i=1:n b(i)=0; for j=1:n b(i)=1/(i+j-1)+b(i); endend修改n分别为5、10、20，运行程序能得到相应A、B、a、b。2、分别求系数矩阵A，B的2-条件数利用自带函数求解

5、n=5时cond(A,2)= 5.4177e+005cond(B,2)= 4.7230e+005n=10时cond(A,2)= 8.6666e+011cond(B,2)= 1.3655e+013n=2=时cond(A,2)= 3.6488e+022cond(B,2)= 1.7243e+018典型病态方程3、利用LU分解法解方程组首先，编辑一个LU分解函数如下functionL,U=Lu(A)% 求解线性方程组的三角分解法% A为方程组的系数矩阵%L和U为分解后的下三角和上三角矩阵n,m=size(A);if n=m error(The rows and columns of matrix A

6、must be equal!); return;end%判断矩阵能否LU分解for ii=1:n for i=1:ii for j=1:ii AA(i,j)=A(i,j); end endif (det(AA)=0) error(The matrix can not be divided by LU!) return;endend%开始计算，先赋初值L=eye(n);U=zeros(n,n);%计算U的第一行，L的第一列for i=1:n U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1);end%计算U的第r行，L的第r列for i=2:n for j=i:n for k

7、=1:i-1 M(k)=L(i,k)*U(k,j); end U(i,j)=A(i,j)-sum(M); end for j=i+1:n for k=1:i-1 M(k)=L(j,k)*U(k,i); end L(j,i)=(A(j,i)-sum(M)/U(i,i); endend然后，编辑一个通过LU分解法解线性方程组的函数如下function L,U,x=Lu_x(A,d)%三角分解法求解线性方程组，LU法解线性方程组Ax=LUx=d%A为方程组的系数矩阵%d为方程组的右端项%L和U为分解后的下三角和上三角矩阵%x为线性方程组的解n,m=size(A);if n=m error(The r

8、ows and columns of matrix A must be equal!); return;end%判断矩阵能否LU分解for ii=1:n for i=1:ii for j=1:ii AA(i,j)=A(i,j); end endif (det(AA)=0) error(The matrix can not be divided by LU!) return;endendL,U=Lu(A); %直接调用自定义函数，首先将矩阵分解，A=LU%设Ly=d由于L是下三角矩阵,所以可求y(i)y(1)=d(1);for i=2:n for j=1:i-1 d(i)=d(i)-L(i,j)

9、*y(j); end y(i)=d(i);end%设Ux=y,由于U是上三角矩阵，所以可求x(i)x(n)=y(n)/U(n,n);for i=(n-1):-1:1 for j=n:-1:i+1 y(i)=y(i)-U(i,j)*x(j); end x(i)=y(i)/U(i,i);end然后，n=5时，调用自定义函数 L,U,x=Lu_x(A,a)解出：x =0.9655 1.2609 0.1625 1.9984 0.6114 L,U,x=Lu_x(B,b)解出：x =0.9976 1.0342 0.8819 1.1477 0.9388第三章 解线性方程组的迭代法第7题有线性方程组，分别写出

10、Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵、收敛性【解】1、Jacobi迭代法Jacobi迭代法的计算公式有迭代矩阵为利用matlab计算pr=max(abs(eig(B)得出迭代矩阵谱半径pr=1.7552 1迭代法不收敛。2、Gauss-Seidel迭代法Gauss-Seidel迭代法的计算公式有迭代矩阵为利用matlab计算pr=max(abs(eig(B)得出迭代矩阵谱半径pr= 0.2746 1迭代法收敛。第9题有方程组，分别写出Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵。用迭代收敛的充要条件给出两种迭代法都收敛的a的取值范围。【解】1、J

11、acobi迭代法Jacobi迭代法的计算公式有迭代矩阵为迭代收敛的充要条件迭代矩阵谱半径小于1即：首先，求出迭代矩阵B的特征值迭代矩阵B的特征方程为得出特征值 代入迭代收敛充要条件得出2、Gauss-Seidel迭代法Jacobi迭代法的计算公式有以上公式矩阵表示为其中可求得迭代矩阵迭代收敛的充要条件迭代矩阵谱半径小于1即：首先，求出迭代矩阵B的特征值迭代矩阵B的特征方程为得出特征值 代入迭代收敛充要条件得出第14题试分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法解线性方程组迭代初始向量取【解】1、Jacobi迭代法（1）编写一个Jacobi迭代法函数，用来输入系数矩阵和

12、右端项，输出解向量，如下functionx,k=Jacobi(A,b,x0,eps,M)%雅可比迭代法求方程组的解%A为方程组的系数矩阵；b为方程组的右端项，列矩阵%x为线性方程组的解，列矩阵；x0为迭代初值，列矩阵%eps为误差限；%M为迭代的最大次数%k当前迭代次数if nargin=3 eps= 1.0e-6; %默认精度 M = 10000; %参数不足时默认后两个条件elseif nargin =4 M = 10000; %参数的默认值elseif nargin=1 error(迭代矩阵谱半径大于1迭代法不收敛); return;endk=0;tol=1;while tol=eps

13、x = B*x0+g; k = k+1; %迭代步数 tol = norm(x-x0);%前后两步迭代结果的误差 x0 = x; if(k=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endend（2）输入系数矩阵A、右端项向量b和初始向量x0A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15b=12;-27;14;-17;12x0=0;0;0;0;0（3）按定义格式输入函数得出结果，如下 x,k=Jacobi(A,b,x0)x = 1.7766 -2.9411 2.8256 -1.

14、3379 0.6040迭代次数k =672、Gauss-Seidel迭代法（1）编写一个Gauss-Seidel迭代法函数，用来输入系数矩阵和右端项，输出解向量，如下functionx,k=GaussSeidel(A,b,x0,eps,M)%高斯赛德尔迭代法求方程组的解（矩阵公式求解）%A为方程组的系数矩阵；b为方程组的右端项%x为线性方程组的解了；x0为迭代初值%eps为误差限；M为迭代的最大次数if nargin=3 eps= 1.0e-6;%默认精度 M = 10000;%参数不足时默认后两个条件elseif nargin =4 M = 10000;%参数的默认值elseif nargi

15、n=1 error(迭代矩阵谱半径大于1迭代法不收敛); return;endk=0;tol=1;while tol=eps x = B*x0+g; k = k+1; %迭代步数 tol = norm(x-x0);%前后两步迭代结果的误差 x0 = x; if(k=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endend（2）输入系数矩阵A、右端项向量b和初始向量x0A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15b=12;-27;14;-17;12x0=0;0;0;0;0（3）按

16、定义格式输入函数得出结果，如下 x,k=GaussSeidel(A,b,x0)x = 1.7070 -2.7326 2.3986 -1.6684 0.1532迭代次数k =38第四章 矩阵特征值与特征向量第2题用幂法求矩阵的按模最大特征值及相应的特征向量，取，要求至少迭代6次。【解】1、编写一个幂法函数functionx,r,k=pmethod(A,x0,eps,M)%幂法求主特征值和特征向量%r是特征值，x为对应特征向量%A目标矩阵%x0为初始向量%eps为误差限%M为迭代的最大次数if nargin=2 eps= 1.0e-6; M = 10000;%参数不足时默认后两个条件elseif

17、nargin =3 M = 10000;%参数的默认值elseif nargin=eps a = max(abs(x); y = x/a; x = A*y; %迭代步数 r = a; tol = abs(r-u);%前后两步迭代结果的误差 u = r; k = k+1; if(k=M) disp(Warning: 迭代次数太多，输出失败！); return; endend2、输入已知参量 A=4 -1 1;16 -2 -2;16 -3 -1 x0=0.5;0.5;13、按定义格式调用函数 x,r,k=pmethod(A,x0,1.0e-6,10)4、输出结果特征向量x = 1.3745 3.6

18、629 3.5804特征值r = 4.0425迭代次数k = 10第4题求矩阵的接近9.6的特征值及相应的特征向量。【解】1、编写一个带原点移位的反幂法函数functionx,r,k=p_method(A,x0,r0,eps,M)%带原点移位的反幂法求按模最小特征值和特征向量%r是按模最小特征值，x为相应特征向量%A目标矩阵；x0为初始向量；r0为某特征值的近似值%eps为误差限；M为迭代的最大次数if nargin=2 r0= 0; %若不指定，则指定其为0 eps= 1.0e-6;%若未指定，则取默认值 M = 10000; %若未指定，则取默认值elseif nargin =3 eps=

19、 1.0e-6;%若未指定，则取默认值 M = 10000; %若未指定，则取默认值elseif nargin =3 M = 10000; %若未指定，则取默认值elseif nargin=eps y = x/a; L,U,x=Lu_x(B,y); %迭代步数 a = max(abs(x); b = a; tol = abs(1/b)-(1/u);%前后两步迭代结果的误差 u = b; k = k+1; if(k=M) disp(Warning: 迭代次数太多，输出失败！); return; endendr = r0 +(1/b);2、输入已知参量 A=1 2 3;2 3 4;3 4 5 r0

20、=9.6 x0=0;0;13、按定义格式调用函数 x,r,k=p_method(A,x0,r0)4、输出结果特征向量x = 22.6520 32.4976 42.1189特征值r = 9.7484迭代次数k = 4第13题已知矩阵试用幂法求按模最大的特征值与特征向量【解】1、输入已知参量 A=190 66 -84 30;66 303 42 -36;336 -168 147 -112;30 -36 28 291 x0=0;0;0;12、按定义格式调用函数幂法求按模最大的特征值与特征向量的函数，参见习题四第2题。 x,r,k=pmethod(A,x0)3、输出结果特征向量x = 1.0e+002

21、* -1.5875 -3.6620 -0.0005 1.7308按模最大特征值r = 3.6057e+002迭代次数k = 103第六章 函数逼近第16题实验数据使用次数x容积y使用次数x容积y2106.4211110.593108.2612110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.30选用双曲线对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。【解】1、编写matlab程序直接调用自带最小二乘法拟合函数编写程序，如下clear,clc;%题目条件x=2 3 5 6

22、7 9 10 11 12 14 16 17 19 20;y=106.42,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;%使用最小二乘法求出1次多项式拟合系数a=polyfit(1./x,1./y,1);%绘制拟合图像xx=0.04:0.01:0.5;yy=a(1)*xx + a(2);plot(1./xx,1./yy,x,y,*);hold on;xx=-0.5:0.01:-0.04;yy=a(1)*xx + a(2);plot(1./xx,1./yy);

23、2、运行程序输出结果使用最小二乘法拟合的曲线方程为下图为绘制出的拟合曲线，并同时将一直点用“*”表示到图中。第七章 数值微分与积分第26题【解】1、编写一个复化Simpson法积分函数function S=FSimpson(f,a,b,eps) %利用复化 Simpson 公式计算被积函数 f(x)在给定区间上的积分值% f:被积函数句柄 % a,b:积分区间的两个端点 % n:子区间个数，n为偶数% S:用复化Simpson法求得的积分值% eps:精度要求if a=b S=0; return;endn=2;h=(b-a)/2;fa=feval(f,a); %计算函数在a点的值 fb=fev

24、al(f,b); %计算函数在b点的值x=a+h; fx=feval(f,x);S=(fa+fb+4*fx)*h/3;S0=S+100;while abs(S-S0)=15*eps S0=S; n=n+2; S=fa+fb; x=a; h=(b-a)/n; %计算步长 for i=1:(n-2)/2 x=x+h; fx=feval(f,x); S=S+4*fx; x=x+h; fx=feval(f,x); S=S+2*fx; end x=x+h; fx=feval(f,x); S=S+4*fx; S=h*S/3;end2、编写matlab程序clear,clc;%题目条件%计算x和y,绘制拟合

25、图像syms t;s=-5;for i=1:101 %取-5到5上的101个点作图 %调用编写函数，使用函数句柄 x(i)= FSimpson(t)(cos(t.2)./2),0,s+(i-1)*0.1,10(-3); y(i)= FSimpson(t)(sin(t.2)./2),0,s+(i-1)*0.1,10(-7);endplot(x,y);2、运行程序输出结果第八章 非线性方程解法题目求下列方程的非零根【解】1、确定根存在的小区间对f（x）求一阶导数得得出单调区间x-771.1-629.8629.8771.1f(x)-+-+-通过上表可以看出，x=629.8时函数可以取到一个极小值，x

26、=771.1时函数可以去到极大值。容易算出，f(630)0。所以，在区间(630,770)之间可以取到一个非零根。2、编写对分法matlab程序计算非零根在matlab中写入题目函数，对分区间来求值clear,clc;%题目条件f=(t)(log(513+0.6651*t)/(513-0.6651*t)-(t/128.52);a=630;%给定区间b=770;fa=feval(f,630);fb=feval(f,770);fx=fa;eps=10(-5);%给定精度N=1000; %限制对分次数n=0;while (fx=0)&(b-a)/2=eps) x=(a+b)/2; fx=feval(

27、f,x); if fa*fxN error(对分次数太多,此区间中不一定有根) return endendxx=x %输出xn %输出迭代次数3、运行程序输出结果xx = 7.5774e+002n = 23此结果的误差限为0.00001第九章 常微分方程数值解法第19题有常微分方程初值问题 分别使用经典RK法和四阶Adams预测-校正算法，求解常微分方程数值解，并与其精确解进行比较，输出结果。其中多步法需要的初值由经典RK法提供。【解】1、编写一个经典RK法求解微分方程数值解的函数function R=Rungkuta4(f,a,b,n,ya) % 功能：用四阶 Runge-Kutta 法求解

28、常微分方程% f:微分方程右端函数句柄 % a,b:自变量取值区间的两个端点 % n:区间等分的个数，ya:函数初值 y(a) % R=x,y:自变量 X 和解 Y 所组成的矩阵 h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; for i=1:n k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+h*k1/2); k3=feval(f,x(i)+h/2,y(i)+h*k2/2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(k1+2*k

29、2+2*k3+k4)*h/6; end R=x,y;2、编写一个四阶Adams预测校正算法的函数function A=CAdams4PC(f,a,b,n,ya)% 功能：用四阶 Adams 预报校正系统求解常微分方程% f:微分方程右端函数句柄 % a,b:自变量取值区间的两个端点 % n:区间等分的个数，ya:函数初值 y(a) % A=x,y:自变量 X 和解 Y 所组成的矩阵 if n4 error(区间个数太小); return; end h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; F=zeros(1,4);

30、 for i=1:n if i4 %用经典RK法，求出y2，y3，y4 k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); elseif i=4 F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); p=feval(f,x(i+1),py); F=F(2) F(3) F(4) p; y(i+1)=y(i)+(h/24)*(F*1,-5,19,9); p=py;c=y(i+1); else F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37