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

数值分析———Matlab上机作业学院：班级：老师：姓名：学号：1第二章解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a，b，c，d输出结果x，均为数组形式functionx=Zhuiganfa(a,b,c,d)%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。%定义三对角矩阵A的各组成单元。方程为Ax=d%b为A的对角线元素(1~n)，a为-1对角线元素(2~n)，c为+1对角线元素(1~n-1)。%A=[2-100%-13-20%0-24-3%00-35]%a=[-1-2-3];c=[-1-2-3];b=[2345];d=[61-21];n=length(b);u(1)=b(1);y(1)=d(1);fori=2:nl(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);fori=(n-1):-1:1%Ux=y,由于U是上三角矩阵，所以可求x(i)x(i)=(y(i)-c(i)*x(i+1))/u(i);end2、输入已知参数a=[2222222];b=[25555555];c=[2222222];d=[220/270000000];23、按定义格式调用函数x=zhuiganfa(a,b,c,d)4、输出结果x=[8.147775166909105-4.0737010928350302.036477565178471-1.0174928201111480.507254485099400-0.2506433926373500.119353996493976-0.047741598597591]第15题【解】1、编写一个程序生成题目条件生成线性方程组Ax=b的系数矩阵A和右端项量b，分别定义矩阵A、B、a、b分别表示系数矩阵，其中1(10.1;,1,2,...,)jijiiaxxiijn或1(,1,2,...,)1ijaijnij分别构成A、B对应右端项量分别a、b。程序如下：clear,clc;n=5;%定义A矩阵A=zeros(n,n);fori=1:nx=1+0.1*i;forj=1:nA(i,j)=x^(j-1);endend%定义B矩阵B=zeros(n,n);fori=1:nforj=1:nB(i,j)=1/(i+j-1);endend%定义a向量，其中Ax=afori=1:nx=1+0.1*i;a(i)=0;forj=1:na(i)=x^(j-1)+a(i);3endend%定义b向量，其中Bx=bfori=1:nb(i)=0;forj=1:nb(i)=1/(i+j-1)+b(i);endend修改n分别为5、10、20，运行程序能得到相应A、B、a、b。2、分别求系数矩阵A，B的2-条件数利用自带函数求解n=5时cond(A,2)=5.361484750384177e+005cond(B,2)=4.766072502417230e+005n=10时cond(A,2)=8.682266932586666e+011cond(B,2)=1.602467527403655e+013n=2=时cond(A,2)=3.420511387496488e+022cond(B,2)=1.845813098127243e+018典型病态方程3、利用LU分解法解方程组首先，编辑一个LU分解函数如下function[L,U]=Lu(A)%求解线性方程组的三角分解法%A为方程组的系数矩阵%L和U为分解后的下三角和上三角矩阵[n,m]=size(A);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判断矩阵能否LU分解forii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)4error('ThematrixcannotbedividedbyLU!')return;endend%开始计算，先赋初值L=eye(n);U=zeros(n,n);%计算U的第一行，L的第一列fori=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);end%计算U的第r行，L的第r列fori=2:nforj=i:nfork=1:i-1M(k)=L(i,k)*U(k,j);endU(i,j)=A(i,j)-sum(M);endforj=i+1:nfork=1:i-1M(k)=L(j,k)*U(k,i);endL(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);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判断矩阵能否LU分解forii=1:nfori=1:iiforj=1:ii5AA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endend[L,U]=Lu(A);%直接调用自定义函数，首先将矩阵分解，A=LU%设Ly=d由于L是下三角矩阵,所以可求y(i)y(1)=d(1);fori=2:nforj=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i);end%设Ux=y,由于U是上三角矩阵，所以可求x(i)x(n)=y(n)/U(n,n);fori=(n-1):-1:1forj=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end然后，n=5时，调用自定义函数[L,U,x]=Lu_x(A,a)解出：x=0.9999999999896551.0000000000326090.9999999999616251.0000000000199840.999999999996114[L,U,x]=Lu_x(B,b)解出：x=0.9999999999999761.0000000000003420.9999999999988191.0000000000014770.9999999999993886第三章解线性方程组的迭代法第7题有线性方程组，分别写出Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵、收敛性123123123104413410811481025xxxxxxxxx【解】1、Jacobi迭代法Jacobi迭代法的计算公式有(1)()()123(1)()()213(1)()()3120.40.41.30.40.81.10.40.82.5kkkkkkkkkxxxxxxxxx迭代矩阵为00.40.40.400.80.40.80B利用matlab计算pr=max(abs(eig(B)))得出迭代矩阵谱半径pr=1.0928203230275521迭代法不收敛。2、Gauss-Seidel迭代法Gauss-Seidel迭代法的计算公式有(1)()()123(1)(1)()213(1)(1)(1)3120.40.41.30.40.81.10.40.82.5kkkkkkkkkxxxxxxxxx迭代矩阵为00.40.400.160.6400.0320.672B利用matlab计算pr=max(abs(eig(B)))得出迭代矩阵谱半径pr=0.6282639865827461迭代法收敛。7第9题有方程组，分别写出Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵。用迭代收敛的充要条件给出两种迭代法都收敛的a的取值范围。112233141001aaxbaxbaxb【解】1、Jacobi迭代法Jacobi迭代法的计算公式有(1)()()1231(1)()212(1)()3134kkkkkkkxaxaxbxaxbxaxb迭代矩阵为040000aaBaa迭代收敛的充要条件迭代矩阵谱半径小于1即：1()max1iinB首先，求出迭代矩阵B的特征值迭代矩阵B的特征方程为3405550aaIBaaaa得出特征值1025a35a代入迭代收敛充要条件得出5555a2、Gauss-Seidel迭代法Jacobi迭代法的计算公式有8(1)()()1231(1)(1)212(1)(1)3134kkkkkkkxaxaxbxaxbxaxb以上公式矩阵表示为(1)(1)()kkkxLxUxg(1)1()1()()kkxILUxILg其中00040000Laa0000000aaU123gbbb可求得迭代矩阵122220()0440aaBILUaaaa迭代收敛的充要条件迭代矩阵谱半径小于1即：1()max1iinB首先，求出迭代矩阵B的特征值迭代矩阵B的特征方程为2222044550aaIBaaaaaa得出特征值1025a35a代入迭代收敛充要条件得出5555a9第14题试分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法解线性方程组1234510123412191232721735143231211743511512xxxxx迭代初始向量取(0)T(0,0,0,0,0)x【解】1、Jacobi迭代法（1）编写一个Jacobi迭代法函数，用来输入系数矩阵和右端项，输出解向量，如下function[x,k]=Jacobi(A,b,x0,eps,M)%雅可比迭代法求方程组的解%A为方程组的系数矩阵；b为方程组的右端项，列矩阵%x为线性方程组的解，列矩阵；x0为迭代初值，列矩阵%eps为误差限；%M为迭代的最大次数%k当前迭代次数ifnargin==3eps=1.0e-6;%默认精度M=10000;%参数不足时默认后两个条件elseifnargin==4M=10000;%参数的默认值elseifnargin3error('参数不足');returnend[n,m]=size(A);nb=length(b);%当方程组行与列的维数不相等时，停止计算，并输出出错信息ifn~=merror('矩阵A行数和列数必须相等!');return;end%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息ifn~=nberror('矩阵A的行数必须和b的长度相等!');return;10endD=zeros(n,n);fori=1:nifA(i,i)==0error('A对角线元素为零！')return;endD(i,i)=A(i,i);%得到矩阵DendB=inv(D)