《微分方程数值解法实验报告》由会员分享，可在线阅读，更多相关《微分方程数值解法实验报告（22页珍藏版）》请在金锄头文库上搜索。
1、微分方程数值解法实验报告： _一：问题描述求解边值问题：其精确解为问题一：取步长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法
2、与PR法比较精度和迭代步数。二实验目的：分别使用五点差分法（Jacobi迭代，Gauss_Seidel迭代，SOR迭代），PR交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。三实验原理：(1) Jacobi迭代法设线性方程组 (1)的系数矩阵A可逆且主对角元素均不为零,令并将A分解成 (2)从而(1)可写成令其中. (3)以为迭代矩阵的迭代法(公式) (4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为 (5)其中为初始向量.(2) Guass-Seidel迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分
3、量,显然在计算第i个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯塞德（Gauss-Seidel）迭代法.把矩阵A分解成 (6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中 (7)以为迭代矩阵构成的迭代法(公式) (8)称为高斯塞德尔迭代法(公式),用 量表示的形式为(3) SOR迭代(4) 交替方向迭代法（PR法）迭代格式为：对于单参数PR法，对于多参数，(5) 共轭梯度法算法步骤如下：预置步任意，计算，并令取：指定算法终止常数，置，进
4、入主步； 主步（）如果，终止算法，输出；否则下行； （）计算：（）计算：（）置，转入（）(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,
5、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)errflag=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迭代法
6、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
7、 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=zero
8、s(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)*
9、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=
10、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(共轭梯度法迭代差
