数值线性代数第二版徐树方高立张平文上机习题第四章实验报告_第1页
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告_第2页
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告_第3页
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告_第4页
数值线性代数第二版徐树方高立张平文上机习题第四章实验报告_第5页
已阅读5页,还剩4页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

第四章上机习题1考虑两点边值问题容易知道它的精确解为为了把微分方程离散化,把[0,1]区间n等分,令h=1/n,得到差分方程简化为从而离散化后得到的线性方程组的系数矩阵为对分别用Jacobi迭代法,G-S迭代法和SOR迭代法求线性方程组的解,要求有4位有效数字,然后比较与精确解得误差。对考虑同样的问题。解(1)给出算法:为解,令,其中,利用Jacobi迭代法,G-S迭代法,SOR迭代法解线性方程组,均可以下步骤求解:step1给定初始向量x0=(0,0,...,0),最大迭代次数N,精度要求c,令k=1step2令x=B*x0+gstep3若||x-x0||2<c,算法停止,输出解和迭代次数k,否则,转step4step4若k>=N,算法停止,迭代失败,否则,令x0=x,转step2在Jacobi迭代法中,B=D-1*(L+U),g=D-1*b在G-S迭代法中,B=D-1*(L+U),g=D-1*b在SOR迭代法中,B=(D-w*L)-1*[(1-w)*D+w*U],g=w*(D-w*L)-1*b另外,在SOR迭代法中,上面算法step1中要给定松弛因子w,其中0<w<2为计算结果,规定w=0.5。给定方程在Ax=b中,矩阵A如题目所示,且为n-1阶矩阵由于在第1、n-1个方程中,由于y(0)、y(1)的存在,使方程右端应减去y(0)、y(1)项,即b(1)=a*h2-c*y(0),b(n-1)=a*h2-(c+h)*y(n),b(i)=a*h2,i=2,...,n-1程序为1Jacobi迭代法编成的函数[x,k]=Jacobi(A,b,c,N)function[x,k]=Jacobi(A,b,c,N)D=diag(diag(A));L=triu(A)-A;U=tril(A)-A;B=D^(-1)*(L+U);g=D^(-1)*b;x0=zeros(length(A),1);x=B*x0+g;k=1;whilenorm(x-x0,2)>=0.00001x0=x;x=B*x0+g;k=k+1;ifk>=Nbreakendendend2G-S迭代法编成的函数[x,k]=GaussSeidel(A,b,c,N)function[x,k]=GaussSeidel(A,b,c,N)U=diag(diag(A))-triu(A);x0=zeros(length(A),1);B=tril(A)^(-1)*U;g=tril(A)^(-1)*b;x=B*x0+g;k=1;whilenorm(x-x0,2)>=0.00001x0=x;x=B*x0+g;k=k+1;ifk>=Nbreakendendend3SOR迭代法编成的函数[x,k]=SOR(A,b,w,c,N)function[x,k]=SOR(A,b,w,c,N)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x0=zeros(length(A),1);B=(D-w*L)^(-1)*((1-w)*D+w*U);g=w*(D-w*L)^(-1)*b;x=B*x0+g;k=1;whilenorm(x-x0,2)>=0.00001x0=x;x=B*x0+g;k=k+1;ifk>=Nbreakendendend4问题1求解ex4_1clear;clc;%c=1;%c=0.1%c=0.01;c=0.0001;a=1/2;n=100;h=1/n;w=1/2;N=1000000;A=-(2*c+h)*eye(n-1);fori=2:n-1wA(i-1,i)=c+h;A(i,i-1)=c;endb=[a*h^2*ones(n-2,1);a*h^2-(c+h)];fori=1:n-1x(i)=i*h;y(i)=((1-a)/(1-exp(-1/c)))*(1-exp(-x(i)/c))+a*x(i);end[y1,n1]=Jacobi(A,b,c,N);[y2,n2]=GaussSeidel(A,b,c,N);[y3,n3]=SOR(A,b,w,c,N);disp(['c=',num2str(c),'时']);disp(['Jacobi迭代与精确解的差为',num2str(norm(y'-y1,inf))]);disp(['迭代次数为',num2str(n1)]);disp(['G-S迭代与精确解的差为',num2str(norm(y'-y2,inf))]);disp(['迭代次数为',num2str(n2)]);disp(['SOR迭代与精确解的差为',num2str(norm(y'-y3,inf))]);disp(['迭代次数为',num2str(n3)]);计算结果为(1)c=1时Jacobi迭代与精确解的差为0.0021999迭代次数为11796G-S迭代与精确解的差为0.0017027迭代次数为6227SOR迭代与精确解的差为0.004511迭代次数为15367(2)c=0.1时Jacobi迭代与精确解的差为0.0094349迭代次数为5353G-S迭代与精确解的差为0.0093007迭代次数为2797SOR迭代与精确解的差为0.010279迭代次数为7300(3)c=0.01时Jacobi迭代与精确解的差为0.066098迭代次数为532G-S迭代与精确解的差为0.066089迭代次数为318SOR迭代与精确解的差为0.06615迭代次数为834(4)c=0.0001时Jacobi迭代与精确解的差为0.0049526迭代次数为116G-S迭代与精确解的差为0.0049507迭代次数为108SOR迭代与精确解的差为0.0049789迭代次数为267结果分析三种迭代法的误差基本相同,且G-S迭代法的收敛速度明显小于Jacobi迭代法,但SOR迭代法收敛速度较慢,原因是收敛因子非最佳。2考虑偏微分方程其中边界条件为u=1.沿x方向和y方向均匀剖分为N等份,令h=1/N,并设应用中心差分离散化后得到差分方程的代数方程组为取g(x,y)和f(x,y)分别为exp(xy)和x+y,用G-S迭代法求解上述方程组,并请列表比较N=20,40,80时收敛所需要的迭代次数和所用的CPU时间。迭代终止条件为||xk+1-xk||2<10-7.解求解过程同问题1给定方程,将组成(n-1)2维列向量,记为则对Ax=b,由于g(x,y)和f(x,y)分别为exp(xy)和x+y,可得A的形式为由于边值,可得,当中的下标存在0,n时,h2(i*h+j*h)应相应的加上1或2(当有两个变量下标均含有0,n时).程序为(1)Jacobi迭代法的程序同问题1(2)问题2求解程序为ex4_2clear;clc;c=10^(-7);n=[20,40,80];N=1000000;form=1:3h=1/n(m);A=zeros((n(m)-1)^2);b=zeros((n(m)-1)^2,1);fori=1:(n(m)-1)^2ifi>1A(i-1,i)=-1;A(i,i-1)=-1;endifi>n(m)-1A(i,i-n(m)+1)=-1;A(i-n(m)+1,i)=-1;endii=ceil(i/(n(m)-1));ifmod(i,n(m)-1)~=0jj=mod(i,n(m)-1);elsejj=n(m)-1;endA(i,i)=4+exp(ii*jj*h^2);b(i)=h^3*(ii+jj);ifii==1||ii==n(m)-1b(i)=b(i)+1;endifjj==1||jj==n(m)-1b(i)=b(i)+1;endend

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论