数值分析matlab上机作业报告.doc

一、给定向量x0，计算初等反射阵Hk。1程序功能： 给定向量x0，计算初等反射阵Hk。2基本原理： 若的分量不全为零，则由确定的镜面反射阵H使得；当时，由有算法：(1)输入x，若x为零向量，则报错(2)将x规范化，如果M0，则报错同时转出停机否则(3)计算，如果，则(4)(5)计算(6)(7)(8)按要求输出,结束3变量说明：x输入的n维向量；nn维向量x的维数；MM是向量x的无穷范数，即x中绝对值最大的一项的绝对值；p Householder初等变换阵的系数；u Householder初等变换阵的向量Us向量x的二范数；x输入的n维向量；nn维向量x的维数；p Householder初等变换阵的系数；u Householder初等变换阵的向量Uk数k，H*x=y，使得y的第k+1项到最后项全为零；4程序代码：(1)function p,u=holder2(x)%HOLDER2 给定向量x0,计算Householder初等变换阵的p,u%程序功能：函数holder2给定向量x0,计算Householder初等变换阵的p,u；%输入：n维向量x；%输出：p,u。p是Householder初等变换阵的系数，% u是Householder初等变换阵的向量U。n=length(x); % 得到n维向量x的维数；p=1;u=0; % 初始化p,u；M=max(abs(x); % 得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；if M=0 % 如果x=0，提示出错,程序终止; disp(Error: M=0); return;else x=x/M; % 规范化end;s=norm(x); % 求x的二范数if x(1)n %如果k值溢出，报错； disp(Error: kn);endH=eye(n); % 初始化H，并使H(1:k,1:k)=I；p,u=holder2(x(k:n); % 得到计算Householde初等变换阵的系数、向量U；H(k:n,k:n)=eye(n-k+1)-pu*u; % 计算H(k:n,k:n)=I-pu*u；5使用示例:情形1:X为零向量 x=0,0,0,0; H=holderk(x,1)Error: M=0H = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1情形2：K值溢出： x=1,2,3,4; H=holderk(x,5)Error: kn情形3：K值为1： x=2,3,4,5; H=holderk(x,1)H = -0.2722 -0.4082 -0.5443 -0.6804 -0.4082 0.8690 -0.1747 -0.2184 -0.5443 -0.1747 0.7671 -0.2911 -0.6804 -0.2184 -0.2911 0.6361检验: det(H)ans = -1.0000 H*xans = -7.3485 0.0000 0.0000 0.0000情形4：(1)K值为3： x=4,3,2,1; H=holderk(x,3)H = 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.8944 -0.4472 0 0 -0.4472 0.8944检验: det(H)ans = -1 H*xans = 4.0000 3.0000 -2.2361 0(2)K值为2: x=4,3,2,1; H=holderk(x,2)H = 1.0000 0 0 0 0 -0.8018 -0.5345 -0.2673 0 -0.5345 0.8414 -0.0793 0 -0.2673 -0.0793 0.9604 det(H)ans = -1.0000 H*xans = 4.0000 -3.7417 0.0000 0二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。1程序功能：给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出AQR2基本原理： 任一实列满秩的mn矩阵A，可以分解成两个矩阵的乘积，即AQR，其中Q是具有法正交列向量的mn矩阵，R是非奇异的n阶上三角阵。算法：(1)输入n阶矩阵A(2)对，求Househoulder初等反射阵的。(3)计算上三角阵R，仍然存储在A(4)计算正交阵Q(5)按要求输出,结束3变量说明：A输入的n阶矩阵，同时用于存储上三角阵R；n矩（方）阵A的阶数；QQ是具有法正交列向量的n阶矩阵；p,u向量A（k:n,k），对应初等反射阵的,uk，jj，ii循环变量；t1计算上三角阵R的系数tj；t2计算正交矩阵Q的系数ti；4程序代码：function Q,A=qrhh(A)%QRHH 用Householder变换法对n阶矩阵A作正交分解A=QR；%程序功能：函数qrhh用Householder变换法对矩阵A作正交分解A=QR；%输入：n阶矩阵A；%输出：Q,A。Q是具有法正交列向量的n阶矩阵，% A（即R）是非奇异的n阶上三角阵，仍用输入的矩阵A存储。%引用函数：% holder2；示例 p,u=holder2(x);n,n=size(A);%求矩（方）阵A的阶数；Q=eye(n);%构造正交矩阵Q(1)=I；for k=1:n-1 p,u=holder2(A(k:n,k);%向量A（k:n,k），对应初等反射阵的,u for jj=k:n%计算上三角阵R（仍存贮于A） t1=dot(u,A(k:n,jj)/p;%利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n%计算正交矩阵Q t2=dot(u,Q(ii,k:n)/p; %利用向量内积求和 Q(ii,k:n)=Q(ii,k:n)-t2*u; endend5使用示例:(1)A为3阶矩阵： A=1 2 3; 2 3 0; 3 4 5A = 1 2 3 2 3 0 3 4 5 q,r=qrhh(A)q = -0.2673 0.8729 0.4082 -0.5345 0.2182 -0.8165 -0.8018 -0.4364 0.4082r = -3.7417 -5.3452 -4.8107 0 0.6547 0.4364 -0.0000 0.0000 3.2660检验： q*rans = 1.0000 2.0000 3.0000 2.0000 3.0000 0.00003.0000 4.0000 5.0000(2)A为4阶矩阵： A=1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0A = 1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0 q,r=qrhh(A)q = -0.2582 0.0597 -0.2660 -0.9268 -0.5164 -0.1045 0.8434 -0.1049 -0.7746 -0.2688 -0.4662 0.3323 -0.2582 0.9556 -0.0222 0.1399r = -3.8730 -6.7132 -6.7132 -6.1968 0 4.4647 6.4805 -1.4783 0 -0.0000 -3.3070 -3.0178 0 0.0000 0 -1.8187检验： q*rans = 1.0000 2.0000 3.0000 4.0000 2.0000 3.0000 -0.0000 1.0000 3.0000 4.0000 5.0000 6.0000 1.0000 6.0000 8.0000 0数值求解正方形域上的Poisson方程边值问题用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组写成矩阵形式Au=f。其中一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。1. 基本原理：Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。假设已经计算出第k步迭代的解（i=1，2，n），要求下一步迭代的解（i=1，2，n）。首先，用Gauss-Seidel迭代格式计算然后引入松弛因子，用松弛因子对和作一个线性组合。，i=1,2,n将二者合并成为一个统一的计算公式：2. 算法(1)Gauss-Seidel迭代法引入松弛因子w： 五点格式即为： (2)计算步骤： 第一步：给松弛因子赋初值w=1.11.8,给场值u和场源b赋初值 第二步：用不同的w进行迭代计算。置error=0;计算在计算机上采用动态计算形式如果|du|error则error=|du|,如果errore,则停机，输出结果u,k.第三步：比较不同的w的迭代次数，用kk存放最小迭代次数，用ww和uu存放相应的w及u。3 程序 u,k=SOR(u,b,w) %（被下面程序调用）%输入场初值u0、场源b及松弛因子w，通过五点差分格式进行迭代运算，%如果第k+1次的迭代结果与第k次的差小于精度，则可以近似认为第k+1次的迭代%结果是精确解，然后返回迭代次数k和迭代解function u,k=SOR(u,b,w) %输出迭代结果u，及迭代次数km=length(u); %m为u的维数h=1/(m-1); %h为步长N=10000;e=0.0000001; %e为精度for k=1:N %k为记录迭代次数 error=0; for j=2:m-1 for jj=2:m-1 sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1); du=w*(h2*b(jj,j)-sum)/4; %计算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if errorabs(du), error=abs(du); end; %统计误差 end end if errork, kk=k;ww=w(i);uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uuendt=toc %统计程序运算时间 4.计算结果： format short n=10; kk,ww,uu=SOR_5dianchafen(n) t = 0.0310kk = 48ww = 1.6000uu = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 contourf (uu, DisplayName, uu); figure(gcf)图一 超松弛二 用块Jacobi迭代法求解线性方程组Au=f。1. 基本原理：对A做自然分解A=D-L-U=D-(L+U)其中D是有A的对角线元素组成的矩阵，L是由A的对角线以下元素组成的矩阵，U是由A得对角线以上元素组成的矩阵。于是将M=D,N=L+U,代入得到Dx=(L+U)x+b任取x的初值进行迭代2. 算法：（1）Gauss-Sediel迭代法原理五点差分格式： 因为A可以写成块状，即： 如果把每一条线上的节点看作一个组，可以把Au=f表示成块状求解：（2）计算步骤： 第一步：给场值u和场源b赋初值，及定义 第二步：用公式，进行迭代计算 第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。3. 程序k,u=kuai_GaussSeidel(n)%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题%输入n,对x、y轴进行n等分；先确定场u的边界及场源b；%用k和u分别存放迭代次数和精确解function k,u=kuai_GaussSeidel(n) %对x、y轴进行n等分h=1/n; %步长u=zeros(n+1,n+1); %对u赋初值u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %计算场源bfor i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endA=zeros(n-1,n-1); %定义矩阵的子块Afor i=1:n-1 if i1, A(i,i-1)=-1; end if i2&jn, u(2:n,j)=pinv(A)*(u(2:n,j-1)+u(2:n,j+1)+b(2:n,j);end end error=max(max(abs(u-ub); %error是前后两次迭代结果的对应元素的最大误差 if error format short n=10; k,u=kuai_GaussSeidel(n)t = 1.3280k = 93u = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 contourf (u, DisplayName, u); figure(gcf)图二 块Gauss-Sediel迭代法三 （预条件）共轭斜量法求解线性方程组Au=f。1.基本原理（1）预条件共轭斜量法原理（2）预优矩阵的选取 2. 算法3.程序%用J-PCG求解正方形域上的Poisson方程边值问题%输入n,对x、y轴进行n等分；先确定场u的边界及场源b；%用k和u分别返回迭代次数和精确解function k,u=J_CG(n)h=1/n; %h步长u=zeros(n+1,n+1); %对u赋初值u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %计算场源bfor i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1); b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endz=zeros(n-1,n-1);r=zeros(n-1,n+1);p=zeros(n-1,n-1);bb=zeros(n-1,n-1);A=zeros(n-1,n-1); %定义矩阵的子块Afor i=1:n-1 if i1, A(i,i-1)=-1; end if i2&j1&i1&in-1,r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i+1)-p(:,i-1);end if i=n-1,r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i-1);end end z(:,1:n-1)=pinv(A)*r(:,1:n-1); % kk=0;mm=0; for i=1:n-1 kk=kk+dot(z(:,i),r(:,i); mm=mm+dot(zz(:,i),rr(:,i); end beita=kk/mm; %beita是 p=z+beita*p; %对p进行更新，确定搜索方向 error=sqrt(norm(u-ub); %error是统计误差 if error format short n=10; k,u=J_CG(n)t = 0.1100k = 37u = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0037 1.0023 1.0000 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0073 1.0045 1.0000 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0107 1.0066 1.0000 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0136 1.0084 1.0000 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0160 1.0100 1.0000 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0174 1.0110 1.0000 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0175 1.0113 1.0000 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0155 1.0103 1.0000 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 contourf (u, DisplayName, u); figure(gcf) 图三 J-PCG四 三种迭代法效率分析 有场的等值图可以看出，三种迭代方法的结果（达到相同的精度）比较一致，但是J-PCG只需要37次（耗时0.1100秒）迭代即可达到比较好的结果；超松弛迭代法需要48次（耗时0.0310秒）迭代即可达到满意的结果；块Gauss-Sediel迭代法需要93次（耗时1.3280）迭代才到达满意的结果。所以对此问题来说，超松弛迭代法和J-PCG法的效率要好于块Gauss-Sediel迭代法。一. 任务： 用MATLAB语言编写连续函数最佳平方逼近的算法程序（函数式M文件）。并用此程序进行数值试验，写出计算实习报告。二. 程序功能要求： 在书中Page355或Page345的程序leastp.m（见附一）的基础上进行修改，使其更加完善。要求算法程序可以适应不同的具体函数，具有一定的通用性。所编程序具有以下功能：1. 用Lengendre多项式做基，并适合于构造任意次数的最佳平方逼近多项式。可利用递推关系 2. 被逼近函数f(x)不用内联函数构造，而改用M文件建立数学函数。这样，此程序可通过修改建立数学函数的M文件以适用不同的被逼近函数（要学会用函数句柄）。3. 要考虑一般的情况。因此，程序中要有变量代换的功能。4. 计算组合系数时，计算函数的积分采用变步长复化梯形求积法（见附三）。5. 程序中应包括帮助文本和必要的注释语句。另外，程序中也要有必要的反馈信息。6. 程序输入：（1）待求的被逼近函数值的数据点(可以是一个数值或向量)（2）区间端点：a,b。7. 程序输出：（1）拟合系数：（2）待求的被逼近函数值1. 基本原理：设内积空间，取线性无关组构造的子空间 span,对被逼近函数f(x)不属于构造逼近函数s(x):s(x)= 求出系数c0,c1, ,cn,使f(x)-s(x)22=min,则函数s(x)就是f(x)H在中的最佳平方逼近元。若取正交函数组对权函数(x)满足=dx=（i,j=0,1,，n）以Legendre正交多项式为函数的最佳平方逼近多项式，若f(x)=xcos(x),x-1，1，用p1(x)=1,p2(x)=x,p3(x)=(3x2-1)/2,p4(x)=(5x3-3x)/2构造函数f(x)的三次最佳平方逼近多项式 :s(x)=c1p1(x)+c2 p2(x)+c3 p3(x)+c4 p4(x)2 程序功能：该程序可对被拟合函数做以Legendre多项式为基的最佳平方拟合，拟合数据可以是一个点或一个向量，并绘出该函数与拟合多项式以及泰勒展式在指定区间上的对比图。3 使用说明：程序用法简单，只需输入拟合数据点x，拟合区间a、b的值，以及c,s=leasp(x,a,b),程序即可自动运行。4 源程序：(1) 主程序M文件%LEASTP.M: 以Legendre多项式为基的最佳平方拟合function c,s=leastp(x,a,b)%x:输入拟合点%func:输入拟合函数%c:输出拟合系数%s:输出拟合值if nargin3, disp(必须输入三个参数); return;endff=func;for i=1:4 pp(i)=2/(2*(i-1)+1);endc(1)=quad(func,a,b)/pp(1)*(2/(b-a);c(2)=quad(fp2,a,b, , ,a,b)/pp(2)*(2/(b-a);c(3)=quad(fp3,a,b, , ,a,b)/pp(3)*(2/(b-a);c(4)=quad(fp4,a,b, , ,a,b)/pp(4)*(2/(b-a);s=c(1)+c(2)*p2(x,a,b)+c(3)*p3(x,a,b)+c(4)*p4(x,a,b);x0=a:0.15:b;s0=c(1)+c(2)*p2(x0,a,b)+c(3)*p