偏微分方程数值解第一次大作业

1、偏微分方程数值解第一次大作业作业一 ：用有限差分法计算：要求：1） L随机选择2） 光滑（随机），且3） 统计L2-误差4） 规模尽可能大根据题意：将线段0L均匀分成M段，则网格步长，节点分布为共M+1个，可得原方程的差分形式：其中，下标k表示在处取值，（）表示在与（与）中点处的取值。化简上式可得：写成矩阵形式：，得，计算时，取，可以解得原方程的解析解（精确解）为由于A为对称三对角阵且主对角占优，可以用直接法求解此方程组。此外，此方程组还可以用经典的Jacobi迭代，共轭梯度（CG）法以及预处理CG法分别来求解，并比较它们的统计误差以及计算时间。但是，注意到当M很大时，由于我们选取的函数形式，

2、有，所以对于Jacobi迭代的系数矩阵有（M越大，越趋近于1），也即此时Jacobi迭代的收敛速度很慢，所以此时没有必要进行Jacobi迭代的计算了。程序中只对直接法、共轭梯度（CG）法和预处理CG法的统计误差和计算时间进行了比较。由于系数矩阵A为对称三对角阵，所以有以下算法：记 直接法：，以上过程其实为消元过程，再用回代过程求出方程组的各变量：，共轭梯度法： 令 对直到收敛完成，结束。特别注意：由于A为对称三对角，为了减少计算过程的乘法次数以及计算机的存储空间，采用两个数组来表示矩阵A，则有其中，。预处理共轭梯度法：采用的是Jacobi迭代预处理，取令 对直到收敛完成，结束注意：仍用CG法中

3、提到的方法简化计算，其实就是解线性方程组，为了提高计算速度采用直接法解此方程组。至此，所有方法的算法都已经给出。在程序的调试过程中发现，当M>6000时，CG法的收敛速度已经很慢了，这是由于此时系数矩阵A的非常趋近于1。因此程序中对于CG算法取的M最大值为5000次。上述，分别为A的最大，最小特征值。结果如下所示： 表一、计算时间（s）M*100012345678910直接法0.0020.0010.0010.0010.0010.0010.0010.0010.0020.002CG法2.3710.6726.2037.2067.12/预处理CG0.0040.0140.0350.0950.206

4、0.3730.7841.1963.8834.593表二、计算误差M*100012345678910直接法5.09E-071.80E-079.80E-086.43E-084.49E-083.54E-082.41E-082.63E-082.17E-081.25E-08CG法5.09E-071.80E-079.80E-086.43E-084.49E-08/预处理CG5.09E-071.80E-079.80E-086.43E-084.50E-083.54E-082.42E-082.64E-082.17E-081.27E-08图1.M与计算误差的双对数由表一可知，用直接法来解此问题是最快捷的，其计算时间

5、在M<6000内基本不随M变化，CG法和预处理CG法随着规模的增大逐渐增大，同时预处理CG法比CG法快了许多。当然，由于这个矩阵为对称三对角矩阵，是一个大型稀疏矩阵，CG算法以及预处理CG算法的优势不能体现出来。但是，如果系数矩阵仅仅是一个大型矩阵，CG算法尤其是预处理CG算法的优势就非常明显了。至于Jacobi迭代方法，本人也进行了尝试，发现其在M=500时计算时间就达到了15s；至于M=1000时，计算运行了5min还是不能达到收敛，而且其计算误差也比前三者大了至少一个量级，因此对于大型系数矩阵，Jacobi方法不具有可行性。对于计算精度，从表一可以看出三种方法的误差基本差不多，同时

6、从图1可以看出当M<6000时，计算误差随着M的增大稳定的减小，但是在M>6000以后减小幅度不再稳定，甚至误差还会增大，这是由于计算机的固有计算精度所致。再增大M，舍入误差的影响将会变的明显。再往后继续增加计算规模，计算误差只能维持10-8左右，这与之前做的第三次作业的结果是一致的，都是与计算机步长有关。作业二：同样考虑二维问题：要求：1） 用Jacobi迭代2） 用CG迭代3） 用预处理CG迭代4） 统计误差和计算时间5） 最大规模不小于此类问题在实际中可以当作变导热系数、有内热源的二维稳态导热问题。为解决此问题，做矩形网格的差分格式。取定x轴和y轴方向的步长和，设的四个邻点为

7、，用表示的中点，见下图。原方程等价于在阴影区域上对原方程进行积分，利用梯形公式有类似的有此外，有利用上面的四个积分近似式并用差商代替其中的偏导数，可以得到原方程的差分表达式将网格的顶点先按j由小到大，j相同的按i由小到大的次序排列，这种排列方式叫“自然顺序排列”，上述差分形式可以得到一个线性方程组，写成矩阵形式，则有，其中且 可记其中，系数矩阵A有以下几个特点：1） A是块三对角矩阵，共有5条对角线上有非零元素；2） A是不可约对角占优的（由于）3） A为稀疏矩阵由于A的上述特点，方程组可以用经典的Jacobi迭代，共轭梯度（CG）法以及预处理CG法分别来求解，并比较它们的统计误差以及计算时间

8、。需要注意的是，由于A不再是对称三对角矩阵了，在编程过程中，迭代过程的算法与第一题的有所区别：CG法及预处理CG法：其中，且时，。类似第一题，由于A为五对角矩阵，为了减少存储空间以及提高计算效率，采用五个数组来表示矩阵A，这样我们节点的最大规模可以达到256*256，但是如果用矩阵来表示A，节点规模最多只能达到64*64。同时，预处理CG法中，解不能使用直接法（LU分解），而用CG法来代替（也可用FFT技巧等，但是这些技巧需要构造系数矩阵A，相比我的算法直接用五个数组代替稀疏矩阵A效果不明显，此外可取的最大规模小了许多）。但是结果表明，在预处理过程中，使用CG法来解线性方程组并不是很理想，原因

9、在下面会分析。Jacobi迭代：易知Jacobi迭代的格式为由于系数矩阵A为块三对角矩阵，上式可以简化为：其中 或者或者结合第一题的算法以及上述改动，可以得到本题的具体算法了。计算时，取，原方程的精确解为其中，。在程序的调试过程中发现，当时Jacobi迭代法的收敛速度已经很慢了，通过计算表明此时系数矩阵A的谱半径。因此程序中对于Jacobi算法取的M最大值为128次，相应的系数矩阵A为阶。当时，CG算法和预处理CG算法所需的收敛时间已经大于2min，收敛速度这么慢的原因是：随着M的增大，系数矩阵A的逐渐趋近于1，导致迭代的收敛变慢；同时，每一步迭代里所需要的乘法次数随着M线性增加。二者的共同作

10、用，导致收敛时间随着M的增加而急剧增加。程序中对于CG算法和预处理CG法取的M最大值为256次，相应的系数矩阵A为阶，已经远大于所要求的的计算规模。上述，分别为A的最大，最小特征值。统计结果如下：表三、计算时间（s）8*2M12345CG法0.0390.2721.92516.504161.322预处理CG0.0380.2371.91116.344150.861Jacobi0.0330.4085.896105.289/表四、计算误差8*2M12345CG法5.98E-052.99E-051.50E-057.48E-063.74E-06预处理CG5.98E-052.99E-051.50E-057.

11、48E-063.74E-06Jacobi5.98E-052.99E-051.49E-057.17E-06/图2. M与计算误差的双对数图由表三可知， CG法和预处理CG法的计算随着规模的增大逐渐增大，但与第一题的结果不同的是预处理CG比CG法快的并不明显，但是当计算的规模逐渐增大时，预处理CG法的优势逐渐体现出来。造成这种差别的原因有二：一、 在预处理过程中，解不再使用直接法（LU分解），而用CG法来代替，由第一题的结果可知CG法的计算耗时远大于直接法，特别是在矩阵A的规模很大时，这种差别更加明显。但是，由于此时矩阵A不是对称三对角，此时比较简单的直接法不再适用，只能退一步采用效率较低的CG算

12、法。二、 系数矩阵A为块对称三对角矩阵，其有更多的非零元素，增加了每一次迭代步中需要进行的乘法运算次数，从而也降低了运算速率。至于Jacobi迭代方法，发现其在M=128时计算时间超过1min；至于M=256时，计算运行了5min还是不能达到收敛，而且其计算误差也比前二者大了许多，因此对于大型的系数矩阵，Jacobi方法不具有可行性。对于计算精度，从表四可以看出CG法和预处理CG法的误差基本差不多，同时从图2可以看出两种方法的计算误差随着M的增大而减小，由于计算时间太长的限制，没有继续增大M以观察其规律，但可以预见是由于计算机的计算精度所限计算误差不可能随着M的增大一直减小下去。再观察Jaco

13、bi迭代，可以发现在时其计算误差随着计算规模的增大而不断减小，而且甚至比CG法和预处理CG法误差还小；但是若再增大计算规模，当时计算误差开始变得很大（结果没有列入表中），结果已经没有可信度了。原因其实很简单，随着计算规模的增大，达到收敛所需的次数急剧增加，而每次迭代的舍入误差不断积累，导致最终的收敛值与真值存在很大的偏差。这从另一方面也说明了对于大型矩阵，Jacobi方法不具有可行性。相反的，CG法和预处理CG法则能够用于求解巨型方程组，但是也有限制条件系数矩阵A必须是对称正定矩阵。另外，在作业过程中，原想取这样一个简单的函数，但是发现M与计算误差的双对数图如下：特别的，当时，计算误差竟然小于

14、；M较小时，误差先随着M不断增大，然后又开始呈下降趋势。这与我们认为的随着节点越多误差越小的事实相矛盾，但其实出现这种现象是合理的，由于u为简单多项式，因此理论上精确解等于差分近似解，误差是计算精度以及舍入误差共同作用的结果。因此，取的函数u为指数函数这种无限光滑的函数形式，而不是简单的多项式形式。附录第一题程序：主程序直接法，返回近似解L=1;format long;M_max=10;err=zeros(3,M_max);t=zeros(3,M_max);k=zeros(2,M_max);for z=1:M_max M=z*1000; h=L/M; a=zeros(M-1,1); b=zer

15、os(M-2,1); c=zeros(M-1,1); for i=1:M-2 a(i)=(exp(i+1/2)*h)+exp(i-1/2)*h)/h2+2*exp(i*h); b(i)=-1/h2*exp(i+1/2)*h); end a(M-1)=(exp(M-1+1/2)*h)+exp(M-1-1/2)*h)/h2+2*exp(M-1)*h); g=2*ones(M-1,1);% 原问题的精确解 x_e1=zeros(M+1,1); c1=(1-exp(1)/(exp(3)-1);c2=(exp(1)-exp(3)/(exp(3)-1); for i=1:M+1 x=1/M*(i-1);

16、x_e1(i)=c1*exp(x)+c2*exp(-2*x)+exp(-x); end x_e=x_e1(2:M);% 直接法求解AU=g t1=clock; x=lu_method(a,b,g,M-1); err(1,z)=norm(x-x_e); t2=clock; t(1,z)=etime(t2,t1);% CG法 if (z<=5) x_cg,k(1,z)=CG_method(a,b,g,M-1); t3=clock; t(2,z)=etime(t3,t2); err(2,z)=norm(x_cg-x_e); else t(2,z)=0; end% 预处理CG法 x_pcg,k(

17、2,z)=PCG_method(a,b,g,M-1); t4=clock; t(3,z)=etime(t4,t3); err(3,z)=norm(x_pcg-x_e);endfunction x=lu_method(a,b,f,n)u=zeros(n-1,1);q=zeros(n,1);u(1)=b(1)/a(1);q(1)=f(1)/a(1);for i=2:n-1 u(i)=b(i)/(a(i)-u(i-1)*b(i-1);endfor i=2:n q(i)=(f(i)-q(i-1)*b(i-1)/(a(i)-u(i-1)*b(i-1);endx=zeros(n,1);x(n)=q(n);

18、for k=1:n-1 i=n-k; x(i)=q(i)-u(i)*x(i+1);endCG法，返回近似解和迭代次数function x1,k=CG_method(a,b,g,m)x0=zeros(m,1);r0=g;p0=r0;k=1;r1=zeros(m,1);while norm(r0)>1e-6 k=k+1; p_A=0; for n=1:m-1 p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1); end p_A=p_A+a(m)*p0(m)*p0(m); ar=(r0'*r0)/p_A; x1=x0+ar*p0; r1(1)=r

19、0(1)-ar*(a(1)*p0(1)+b(1)*p0(2); for i=2:m-1 r1(i)=r0(i)-ar*(b(i-1)*p0(i-1)+a(i)*p0(i)+b(i)*p0(i+1); end r1(m)=r0(m)-ar*(a(m)*p0(m)+b(m-1)*p0(m-1); p0=r1+(r1'*r1)/(r0'*r0)*p0; r0=r1; x0=x1;endend预处理CG法，返回近似解和迭代次数function x1,k=PCG_method(a,b,g,m)x0=zeros(m,1);r0=g;z1=lu_method(a,b,r0,m);p0=z1;

20、k=1;r1=zeros(m,1);while norm(r0)>1e-8 k=k+1; p_A=0; for n=1:m-1 p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1); end p_A=p_A+a(m)*p0(m)*p0(m); ar=(r0'*z1)/p_A; x1=x0+ar*p0; r1(1)=r0(1)-ar*(a(1)*p0(1)+b(1)*p0(2); for i=2:m-2 r1(i)=r0(i)-ar*(b(i-1)*p0(i-1)+a(i)*p0(i)+b(i)*p0(i+1); end r1(m)=r0(m)

21、-ar*(a(m)*p0(m)+b(m-1)*p0(m-1); z1=lu_method(a,b,r1,m); p0=z1+(r1'*z1)/(r0'*z1)*p0; r0=r1; x0=x1;endend第二题程序：主程序CG法，返回近似解和迭代次数L1=1;L2=1;format long;M_max=5;err=zeros(3,M_max);t=zeros(3,M_max);k=zeros(2,M_max);c1=1/(exp(1)+1);c2=exp(1)/(exp(1)+1);for z=1:M_max M=8*2z; h1=L1/M; h2=L2/M; a=zero

22、s(M-1)2,1); b=zeros(M-1)2-1,1); c=zeros(M-1)*(M-2),1); f=zeros(M-1)2,1); for j=1:M-1 for i=1:M-1 a(M-1)*(j-1)+i)=(exp(-(i-1/2)*h1)+exp(-(i+1/2)*h1)/h1/h1. +2*exp(-i*h1)/h2/h2+exp(-i*h1); f(M-1)*(j-1)+i)=exp(-i*h1)-c1-2*c2*c2*exp(-2*i*h1-j*h2). +c2*exp(-2*i*h1)-2*c1*c2*exp(j*h2-2*i*h1); if i<M-1 b

23、(M-1)*(j-1)+i)=-exp(-(i+1/2)*h1)/h1/h1; end if j>=2 c(M-1)*(j-2)+i)=-exp(-i*h1)/h2/h2; end end end% 精确解 x_el=zeros(M-1)2,1); for j=1:(M-1) for i=1:M-1 x_el(M-1)*(j-1)+i)=(c1*exp(i*h1)+c2*exp(-i*h1)-1)*. (c1*exp(j*h2)+c2*exp(-j*h2)-1); end end X=zeros(M-1,M-1); for j=1:M-1 for i=1:M-1 X(i,j)=x_el(

24、M-1)*(j-1)+i); end end% CG t1=clock; x_cg,k(1,z)=CG_method(a,b,c,f,M-1); X_cg=zeros(M-1,M-1); for j=1:M-1 for i=1:M-1 X_cg(i,j)=x_cg(M-1)*(j-1)+i); end end t2=clock; t(1,z)=etime(t2,t1); err(1,z)=norm(X_cg-X);% PCG x_pcg,k(2,z)=PCG_method(a,b,c,f,M-1); X_pcg=zeros(M-1,M-1); for j=1:M-1 for i=1:M-1 X

25、_pcg(i,j)=x_pcg(M-1)*(j-1)+i); end end t3=clock; t(2,z)=etime(t3,t2); err(2,z)=norm(X_pcg-X);% Jacobi if M<=200 x_J,k(3,z)=Jacobi(a,b,c,f,M-1); X_J=zeros(M-1,M-1); for j=1:M-1 for i=1:M-1 X_J(i,j)=x_J(M-1)*(j-1)+i); end end t4=clock; t(3,z)=etime(t4,t3); err(3,z)=norm(X_J-X); endendfunction y1,k=

26、CG_method(a,b,c,g,M)m=M2;y0=zeros(m,1);r0=g;p0=r0;k=1;y1=zeros(m,1);r1=zeros(m,1);while norm(r0)>1e-6 k=k+1; p_A=0; for n=1:m-1 if n<=m-M p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1). +2*c(n)*p0(n)*p0(n+M); else p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1); end end p_A=p_A+a(m)*p0(m)*p0(m);

27、ar=(r0'*r0)/p_A; y1=y0+ar*p0; r1(1)=r0(1)-ar*(a(1)*p0(1)+b(1)*p0(2)+c(1)*p0(1+M); for i=2:m-1 if i>M&&i<=m-M r1(i)=r0(i)-ar*(c(i-M)*p0(i-M)+b(i-1)*p0(i-1). +a(i)*p0(i)+b(i)*p0(i+1)+c(i)*p0(i+M); elseif i<=M r1(i)=r0(i)-ar*(b(i-1)*p0(i-1)+a(i)*p0(i). +b(i)*p0(i+1)+c(i)*p0(i+M); e

28、lse r1(i)=r0(i)-ar*(c(i-M)*p0(i-M)+b(i-1)*p0(i-1). +a(i)*p0(i)+b(i)*p0(i+1); end end r1(m)=r0(m)-ar*(a(m)*p0(m)+b(m-1)*p0(m-1)+c(m-M)*p0(m-M); p0=r1+(r1'*r1)/(r0'*r0)*p0; r0=r1; y0=y1;end预处理CG法，返回近似解和迭代次数function x1,k=PCG_method(a,b,c,g,M)m=M2;x0=zeros(m,1);r0=g;z1,=CG_method(a,b,c,r0,M);p0=

29、z1;k=1;r1=zeros(m,1);while norm(r0)>1e-4 k=k+1; p_A=0; for n=1:m-1 if n<=m-M p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1). +2*c(n)*p0(n)*p0(n+M); else p_A=p_A+a(n)*p0(n)*p0(n)+2*b(n)*p0(n)*p0(n+1); end end p_A=p_A+a(m)*p0(m)*p0(m); ar=(r0'*z1)/p_A; x1=x0+ar*p0; r1(1)=r0(1)-ar*(a(1)*p0(1)+b(1)*p0(2)+c(1)*p0(1+M); fo