微分方程数值解法实验报告.docx

微分方程数值解法实验报告姓名： 班级：学号：一：问题描述求解边值问题：其精确解为问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi迭代法，Gauss_Seidel迭代法，SOR迭代法（w=1.45）。求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解与精确解的精度。问题二：取步长h=k=1/64,1/128,作五点差分格式，用单参数和双参数PR法解差分方程，近似到小数点后四位。与SOR法比较精度和迭代步数。问题三：取步长h=k=1/64,1/128,作五点差分格式，用共轭梯度法和预处理共轭梯度法解差分方程，近似到小数点后四位。与SOR法与PR法比较精度和迭代步数。二实验目的：分别使用五点差分法（Jacobi迭代，Gauss_Seidel迭代，SOR迭代），PR交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。三实验原理：(1) Jacobi迭代法设线性方程组 (1)的系数矩阵A可逆且主对角元素均不为零,令 并将A分解成 (2)从而(1)可写成 令 其中. (3)以为迭代矩阵的迭代法(公式) (4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为 (5)其中为初始向量.(2) Guass-Seidel迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯塞德（Gauss-Seidel）迭代法.把矩阵A分解成 (6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中 (7)以为迭代矩阵构成的迭代法(公式) (8)称为高斯塞德尔迭代法(公式),用 量表示的形式为 (3) SOR迭代 (4) 交替方向迭代法（PR法）迭代格式为：对于单参数PR法，对于多参数，(5) 共轭梯度法 算法步骤如下：预置步任意，计算，并令取：指定算法终止常数，置，进入主步； 主步（）如果，终止算法，输出；否则下行； （）计算：（）计算： （）置，转入（） (6) 预共轭梯度法预置步任意，计算，并令取：指定算法终止常数，置，进入主步； 主步（）计算：， （）如果，转入（3）否则，终止算法，输出计算结果 （）计算： （）置，转入（1） 注：在算法主步中，引入变量，及，可以简化计算。四程序设计(MATLAB实现)Jacobi迭代法functionu_1,m_1=Jacobi_Solve(A,b,n,err)D=diag(A);D=diag(D);L=-tril(A,-1);R=-triu(A,1);B=D(L+R);g=(Db);m_1=0;u_1_0=zeros(n-1,n-1);%初始迭代值u_1_0=u_1_0(:);flag=1;while flag u_1=B*u_1_0+g; if norm(u_1-u_1_0,inf)err flag=0; end u_1_0=u_1; m_1=m_1+1;%迭代次数值enduu=zeros(n+1);for mm=1:n-1 for nn=1:n-1 uu(nn+1,mm+1)=u_1(mm-1)*(n-1)+nn,1); endend%Jacobi迭代差分解图像x=0:1/n:1;y=0:1/n:1;figure(2)mesh(x,y,uu);title(Jacobi迭代差分解图像)Gauss_Seidel迭代法functionu_2,m_2=Gauss_Seidel_Solve(A,b,n,err)if nargin=2 err=1e-6;endD=diag(A);D=diag(D);L=-tril(A,-1);U=-triu(A,1);R=(D-L)U;g=(D-L)b;m_2=0;u_2_0=zeros(n-1,n-1);%初始迭代值u_2_0=u_2_0(:);flag=1;while flag u_2=R*u_2_0+g; if norm(u_2-u_2_0,inf)err flag=0; end u_2_0=u_2; m_2=m_2+1;%迭代次数值enduu=zeros(n+1);for mm=1:n-1 for nn=1:n-1 uu(nn+1,mm+1)=u_2(mm-1)*(n-1)+nn,1); endend%Gauss-Seidel迭代差分解图像x=0:1/n:1;y=0:1/n:1;figure(3)mesh(x,y,uu); title(Gauss-Seidel迭代差分解图像)SOR迭代法functionu_3,m_3=SOR_Solve(A,b,err,w,n)D=diag(A);D=diag(D);L=-tril(A,-1);R2=-triu(A,1);B1=(D-w*L)(1-w)*D+w*R2);g=w*(D-w*L)b;m_3=0;u_3_0=zeros(size(b,1),1);%初始迭代值flag=1;while flag u_3=B1*u_3_0+g; if norm(u_3-u_3_0,inf)=0.01 count=count+1; u=u2; u1=inv(I+tao_opt*L_1)*(I-tao_opt*L_2)*u+tao_opt*b); u2=inv(I+tao_opt*L_2)*(I-tao_opt*L_1)*u1+tao_opt*b);endu2z=zeros(N-1)2,1);for j=1:N-1 for i=1:N-1 z(j-1)*(N-1)+i)=exp(pi*(i/N+j/N)*(sin(pi*i/N)*sin(pi*j/N); endendzcountuu=zeros(N+1);for mm=1:N-1 for nn=1:N-1 uu(nn+1,mm+1)=u2(mm-1)*(N-1)+nn,1); endend%画图x=0:1/N:1;y=0:1/N:1;figure(1)mesh(x,y,uu);title(单参数PR迭代差分解图像（c1=c2=0.5*(1/sin(pi/16) %画图数值解共轭梯度法functionu_5,m_5=gongetidufa(A,b,n,err)m_5=0;u_5_0=zeros(size(b,1),1); p0=b-A*u_5_0;r0=p0;flag=1;while flag a0=(r0)*r0)/(p0)*(A*p0); u_5=u_5_0+a0*p0; if norm(u_5-u_5_0,inf)err flag=0; end u_5_0=u_5; r1=r0-a0*A*p0; b0=(r1)*r1)/(r0)*r0); p0=r1+b0*p0; r0=r1; m_5=m_5+1;enduu=zeros(n+1);for mm=1:n-1 for nn=1:n-1 uu(nn+1,mm+1)=u_5(mm-1)*(n-1)+nn,1); endend%共轭梯度法迭代差分解图像x=0:1/n:1;y=0:1/n:1;figure(5)mesh(x,y,uu);title(共轭梯度法迭代差分解图像) 预处理共轭梯度法functionu_6,m_6=yuchuligongetidufa(A,b,n,err)m_6=0;u_6_0=zeros(size(b,1),1); % D=diag(A);% D=diag(D);% L=tril(A,-1);% L=L+DD;% B=L*D*(L);B=0.7*eye(n-1)2)p0=b-A*u_6_0;r0=p0;s0=Br0;flag=1;while flag a0=(r0)*s0)/(p0)*(A*p0); u_6=u_6_0+a0*p0; if norm(u_6-u_6_0,inf)err flag=0; end u_6_0=u_6; r1=r0-a0*A*p0; s1=Br1; b0=(r1)*s1)/(r0)*s0); p0=s1+b0*p0; r0=r1; s0=s1; m_6=m_6+1;enduu=zeros(n+1);for mm=1:n-1 for nn=1:n-1 uu(nn+1,mm+1)=u_6(mm-1)*(n-1)+nn,1); endend%预处理共轭梯度法迭代差分解图像x=0:1/n:1;y=0:1/n:1;figure(6)mesh(x,y,uu);title(预处理共轭梯度法迭代差分解图像) 五实验数据分析1. 实验数据（部分）u精确值u_1u_2u_3u_4u_5u_60.05637 0.04884 0.04891 0.03396 0.04840 0.04901 0.04905 0.13455 0.11953 0.11967 0.08309 0.11869 0.11991 0.11991 0.23772 0.21531 0.21552 0.14965 0.21408 0.21585 0.21586 0.36820 0.33863 0.33891 0.23534 0.33708 0.33930 0.33935 0.52689 0.49057 0.49091 0.34090 0.48875 0.49139 0.49143 0.71247 0.67004 0.67043 0.46557 0.66802 0.67095 0.67094 0.92044 0.87287 0.87330 0.60645 0.87074 0.87383 0.87381 1.14208 1.09068 1.09112 0.75772 1.08851 1.09169 1.09163 1.36315 1.30963 1.31008 0.90978 1.30751 1.31062 1.31056 1.56264 1.50911 1.50954 1.04829 1.50711 1.51002 1.51000 1.71145 1.66037 1.66076 1.15331 1.65858 1.66118 1.66121 1.77124 1.72534 1.72568 1.19840 1.72384 1.72602 1.72607 1.69357 1.65570 1.65597 1.14999 1.65453 1.65624 1.65629 1.41964 1.39248 1.39267 0.96714 1.39166 1.39286 1.39289 0.88074 0.86649 0.86659 0.60181 0.86610 0.86670 0.86671 0.13455 0.11953 0.11967 0.08309 0.11869 0.11991 0.11991 0.32120 0.29120 0.29149 0.20241 0.28962 0.29194 0.29193 0.56747 0.52270 0.52313 0.36326 0.52036 0.52377 0.52378 0.87895 0.81982 0.82038 0.56969 0.81686 0.82120 0.82120 1.25777 1.18506 1.18574 0.82341 1.18158 1.18669 1.18668 1.70076 1.61573 1.61651 1.12257 1.61188 1.61752 1.61749 2.19723 2.10177 2.10261 1.46014 2.09769 2.10368 2.10361 2.72631 2.62302 2.62390 1.82216 2.61887 2.62499 2.62489 3.25404 3.14632 3.14721 2.18557 3.14227 3.14825 3.14817 3.73024 3.62234 3.62319 2.51612 3.61854 3.62413 3.62411 4.08547 3.98234 3.98313 2.76608 3.97894 3.98391 3.98397 4.22818 4.13538 4.13606 2.87229 4.13251 4.13672 4.13678 4.04279 3.96609 3.96663 2.75463 3.96386 3.96714 3.96720 3.38888 3.33378 3.33416 2.31541 3.33224 3.33453 3.33454 2.10245 2.07352 2.07372 1.44009 2.07277 2.07391 2.07391 0.23772 0.21531 0.21552 0.14965 0.21408 0.21585 0.21586 0.56747 0.52270 0.52313 0.36326 0.52036 0.52377 0.52378 1.00258 0.93567 0.93631 0.65019 0.93223 0.93721 0.93721 1.55288 1.46439 1.46523 1.01750 1.46004 1.46636 1.46639 2.22215 2.11316 2.11418 1.46817 2.10805 2.11550 2.11550 3.00480 2.87712 2.87828 1.99880 2.87146 2.87971 2.87969 3.88193 3.73830 3.73956 2.59692 3.73231 3.74106 3.74100 4.81668 4.66094 4.66226 3.23770 4.65485 4.66379 4.66368 5.74904 5.58628 5.5876