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

微分方程数值解法实验报告 姓名： 班级： 学号： 一：问题描述 求解边值问题： 其精确解为 问题一：取步长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迭代法 function[u_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=(D\b); 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;%迭代次数值 end uu=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); end end %Jacobi迭代差分解图像 x=[0:1/n:1]; y=[0:1/n:1]; figure(2) mesh(x,y,uu);title(Jacobi迭代差分解图像) Gauss_Seidel迭代法 function[u_2,m_2]=Gauss_Seidel_Solve(