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

数值分析 ———Matlab上机作业 学院： 班级： 老师： 姓名： 学号： 第二章 解线性方程组的直接解法 第14    【解】 1、编写一个追赶法的函数 输入a，b，c，d输出结果x，均为数组形式 function x=Zhuiganfa(a,b,c,d) %首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。 %定义三对角矩阵A的各组成单元。方程为Ax=d %b为A的对角线元素(1~n)，a为-1对角线元素(2~n)，c为+1对角线元素(1~n-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];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); end x(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); end 2、输入已知参数 >>a=[2 2 2 2 2 2 2]; >>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.  -4.  2.  -1.  0.  -0.  0.  -0.] 第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*i; for j=1:n A(i,j)=x^(j-1); end end %定义B矩阵 B=zeros(n,n); for i=1:n for j=1:n B(i,j)=1/(i+j-1); end end %定义a向量，其中Ax=a for i=1:n x=1+0.1*i; a(i)=0; for j=1:n a(i)=x^(j-1)+a(i); end end      %定义b向量，其中Bx=b for i=1:n b(i)=0; for j=1:n b(i)=1/(i+j-1)+b(i); end end 修改n分别为5、10、20，运行程序能得到相应A、B、a、b。 2、分别求系数矩阵A，B的2-条件数 利用自带函数求解 n=5时 cond(A,2)= 5.e+005 cond(B,2)= 4.e+005 n=10时 cond(A,2)= 8.e+011 cond(B,2)= 1.e+013 n=2=时 cond(A,2)= 3.e+022 cond(B,2)= 1.e+018 典型病态方程 3、利用LU分解法解方程组 首先，编辑一个LU分解函数如下 function[L,U]=Lu(A) % 求解线性方程组的三角分解法 % A为方程组的系数矩阵 %L和U为分解后的下三角和上三角矩阵 [n,m]=size(A); if n~=m error('The rows 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 end if (det(AA)==0) error('The matrix can not be divided by LU!') return; end end %开始计算，先赋初值 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=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); end end 然后，编辑一个通过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 rows 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 end if (det(AA)==0) error('The matrix can not be divided by LU!') return; end end [L,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)*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.  1.  0.  1.  0. >> [L,U,x]=Lu_x(B,b) 解出： x =0.  1.  0.  1.  0. 第三章 解线性方程组的迭代法 第7题 有线性方程组，分别写出Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵、收敛性 【解】 1、Jacobi迭代法 Jacobi迭代法的计算公式有 迭代矩阵为 利用matlab计算 >>pr=max(abs(eig(B))) 得出迭代矩阵谱半径pr=1. > 1迭代法不收敛。 2、Gauss-Seidel迭代法 Gauss-Seidel迭代法的计算公式有 迭代矩阵为 利用matlab计算 >>pr=max(abs(eig(B))) 得出迭代矩阵谱半径pr= 0. < 1迭代法收敛。 第9题 有方程组，分别写出Jacobi和Gauss-Seidel迭代法的计算公式、迭代矩阵。用迭代收敛的充要条件给出两种迭代法都收敛的a的取值范围。 【解】 1、Jacobi迭代法 Jacobi迭代法的计算公式有 迭代矩阵为 迭代收敛的充要条件迭代矩阵谱半径小于1即： 首先，求出迭代矩阵B的特征值 迭代矩阵B的特征方程为 得出特征值 代入迭代收敛充要条件得出 2、Gauss-Seidel迭代法 Jacobi迭代法的计算公式有 以上公式矩阵表示为 其中 可求得迭代矩阵 迭代收敛的充要条件迭代矩阵谱半径小于1即： 首先，求出迭代矩阵B的特征值 迭代矩阵B的特征方程为 得出特征值 代入迭代收敛充要条件得出 第14题 试分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法解线性方程组 迭代初始向量取 【解】 1、Jacobi迭代法 （1）编写一个Jacobi迭代法函数，用来输入系数矩阵和右端项，输出解向量，如下 function[x,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<3 error('参数不足'); return end [n,m]=size(A); nb=length(b); %当方程组行与列的维数不相等时，停止计算，并输出出错信息 if n~=m error('矩阵A行数和列数必须相等!'); return; end %当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if n~=nb error('矩阵A的行数必须和b的长度相等!'); return; end D =zeros(n,n); for i=1:n if A(i,i)==0 error('A对角线元素为零！') return; end D(i,i)=A(i,i);  %得到矩阵D end B=inv(D)*(D-A);      %B为迭代矩阵 g=inv(D)*b;          %g为右端项 pr=max(abs(eig(B))); %求迭代矩阵谱半径 if pr>=1 error('迭代矩阵谱半径大于1迭代法不收敛'); return; end k=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; end end （2）输入系数矩阵A、右端项向量b和初始向量x0 >>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] （3）按定义格式输入函数得出结果，如下 >> [x,k]=Jacobi(A,b,x0) x = 1. -2. 2. -1. 0. 迭代次数k =67 2、Gauss-Seidel迭代法 （1）编写一个Gauss-Seidel迭代法函数，用来输入系数矩阵和右端项，输出解向量，如下 function[x,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 nargin<3 error('参数不足'); return end [n,m]=size(A); nb=length(b); %当方程组行与列的维数不相等时，停止计算，并输出出错信息 if n~=m error('矩阵A行数和列数必须相等!'); return; end %当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if n~=nb error('矩阵A的行数必须和b的长度相等!'); return; end L =zeros(n,n); U =zeros(n,n); D =zeros(n,n); for i=2:n for j=1:i-1 L(i,j)=-A(i,j); end end for i=1:n-1 for j=i+1:n U(i,j)=-A(i,j); end end for i=1:n D(i,i)=A(i,i); end B=inv(D-L)*U;        %B为迭代矩阵 g=inv(D-L)*b;        %g为右端项 pr=max(abs(eig(B))); %求迭代矩阵谱半径 if pr>=1 error('迭代矩阵谱半径大于1迭代法不收敛'); return; end k=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; end end （2）输入系数矩阵A、右端项向量b和初始向量x0 >>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] （3）按定义格式输入函数得出结果，如下 >> [x,k]=GaussSeidel(A,b,x0) x = 1. -2. 2. -1. 0. 迭代次数k =38 第四章 矩阵特征值与特征向量 第2题 用幂法求矩阵 的按模最大特征值及相应的特征向量，取 ，要求至少迭代6次。 【解】 1、编写一个幂法函数 function[x,r,k]=pmethod(A,x0,eps,M) %幂法求主特征值和特征向量 %r是特征值，x为对应特征向量 %A目标矩阵 %x0为初始向量 %eps为误差限 %M为迭代的最大次数 if nargin==2 eps= 1.0e-6; M = 10000;%参数不足时默认后两个条件 elseif nargin ==3 M = 10000;%参数的默认值 elseif nargin<2 error('参数不足'); return end [n,m]=size(A); %当方程组行与列的维数不相等时，停止计算，并输出出错信息 if n~=m error('矩阵A行数和列数必须相等!'); return; end k=0; u=0; tol=1; x=x0; while tol>=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; end end 2、输入已知参量 >> A=[4 -1 1;16 -2 -2;16 -3 -1] >> x0=[0.5;0.5;1] 3、按定义格式调用函数 >> [x,r,k]=pmethod(A,x0,1.0e-6,10) 4、输出结果 特征向量 x = 1. 3. 3. 特征值 r = 4. 迭代次数 k = 10 第4题 求矩阵 的接近9.6的特征值及相应的特征向量。 【解】 1、编写一个带原点移位的反幂法函数 function[x,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= 1.0e-6;%若未指定，则取默认值 M = 10000;  %若未指定，则取默认值 elseif nargin ==3 M = 10000;  %若未指定，则取默认值 elseif nargin<2 error('参数不足'); return end [n,m]=size(A); %当方程组行与列的维数不相等时，停止计算，并输出出错信息 if n~=m error('矩阵A行数和列数必须相等!'); return; end I=eye(n); B=zeros(n,n); B=A-r0*I; k=0; u=1; tol=1; x=x0; a = max(abs(x)); while tol>=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; end end r = r0 +(1/b); 2、输入已知参量 >> A=[1 2 3;2 3 4;3 4 5] >> r0=9.6 >> x0=[0;0;1] 3、按定义格式调用函数 >> [x,r,k]=p_method(A,x0,r0) 4、输出结果 特征向量 x = 22. 32. 42. 特征值 r = 9. 迭代次数 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;1] 2、按定义格式调用函数 幂法求按模最大的特征值与特征向量的函数，参见习题四第2题。 >> [x,r,k]=pmethod(A,x0) 3、输出结果 特征向量 x = 1.0e+002 * -1. -3. -0. 1. 按模最大特征值 r = 3.e+002 迭代次数 k = 103 第六章 函数逼近 第16题 实验数据 使用次数x 容积y 使用次数x 容积y 2 106.42 11 110.59 3 108.26 12 110.60 5 109.58 14 110.72 6 109.50 16 110.90 7 109.86 17 110.76 9 110.00 19 111.10 10 109.93 20 111.30         选用双曲线 对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。 【解】 1、编写matlab程序 直接调用自带最小二乘法拟合函数编写程序，如下 clear,clc; %题目条件 x=[2 3 5 6 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); 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; end n=2; h=(b-a)/2; fa=feval(f,a); %计算函数在a点的值 fb=feval(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; end 2、编写matlab程序 clear,clc; %题目条件 %计算x和y,绘制拟合图像 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)); end plot(x,y); 2、运行程序输出结果 第八章 非线性方程解法 题目 求下列方程的非零根 【解】 1、确定根存在的小区间 对f（x）求一阶导数得 得出单调区间 x   -771.1 -629.8 629.8 771.1   f’(x) - + - + -                       通过上表可以看出，x=629.8时函数可以取到一个极小值，x=771.1时函数可以去到极大值。容易算出，f(630)<0且f(770)>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(f,x); if fa*fx<0 b=x; fb=fx; else a=x; fa=fx; end n=n+1; if n>N error('对分次数太多,此区间中不一定有根') return end end xx=x %输出x n %输出迭代次数 3、运行程序输出结果 xx = 7.e+002 n = 23 此结果的误差限为0.00001 第九章 常微分方程数值解法 第19题 有常微分方程初值问题 分别使用经典RK法和四阶Adams预测-校正算法，求解常微分方程数值解，并与其精确解 进行比较，输出结果。其中多步法需要的初值由经典RK法提供。 【解】 1、编写一个经典RK法求解微分方程数值解的函数 function R=Rungkuta4(f,a,b,n,ya) % 功能：用四阶 Runge-Kutta 法求解常微分方程 % 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*k2+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 n<4 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); for i=1:n if i<4  %用经典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,-59,55]'); my=py-251*(p-c)/270; m=feval(f,x(i+1),my); F=[F(2) F(3) F(4) m]; cy=y(i)+(h/24)*(F*[1,-5,19,9]'); y(i+1)=cy+19*(py-cy)/270; p=py;c=cy; end end A=[x',y']; 3、编写matlab程序按题目条件计算 选取区间个数n=10 clear,clc; %题目条件 f=@(x,y)(-y+2*cos(x)); %使用经典RK法，求数值解 n=10; R=Rungkuta4(f,0,pi,n,1); %使用四阶Adams预测校正算法，求数值解 A=CAdams4PC(f,0,pi,n,1); %求出精确解 f0=@(x)(cos(x)+sin(x)); R0=zeros(1,n+1); R0=feval(f0,R(:,1)); %输出自变量，数值解，精确值，整体阶段误差 RK=[R(:,1),R(:,2),R0,R0-R(:,2)]  %经典RK法结果输出 AD=[A(:,1),A(:,2),R0,R0-A(:,2)]  %四阶Adans预测校正法结果输出 4、运行程序输出结果 RK = 0 1. 1. 0 0. 1. 1. 0. 0.628318********* 1. 1. 0. 0. 1. 1. 0. 1. 1. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 2. 0. 0. 0. 2. -0. -0. 0. 2. -0. -0. 0. 3. -0. -1. -0.         AD = 0 1. 1. 0 0. 1. 1. 0. 0.628318********* 1. 1. 0. 0. 1. 1. 0. 1. 1. 1. 0. 1. 0. 1. 0.         1. 0. 0. 0. 2. 0. 0. 0. 2. -0. -0. 0. 2. -0. -0. 0. 3. -1. -1. 0.