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

微分方程数值解法实验报告姓名：班级：学号：一：问题描述求解边值问题：()2(sincoscossin(0,1)(0,1)0,(,)xyuexyxyGuxyG（x,y)其精确解为)sin()sin(),()(yxeyxuyx问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi迭代法，Gauss_Seidel迭代法，SOR迭代法（w=1.45）。求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解)128/1,64/1)(,(hyxuh与精确解的精度。问题二：取步长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可逆且主对角元素均不为零,令bAxnna,...,a,a2211并将A分解成(2)从而(1)可写成令其中.(3)以为迭代矩阵的迭代法(公式)(4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5)其中为初始向量.(2)Guass-Seidel迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel）迭代法.把矩阵A分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中(7)以为迭代矩阵构成的迭代法(公式)(8)nna,...,a,adiagD2211DDAAbxADDx11fxBxbDf,ADIB11111B111fxBxkk,...,,k,n,...,ixabaxnijj)k(jjiiii)k(i21021111Tnx,...x,xx002010kx1kx1kix1111kikx,...,x1k1kx1kjxULDAnna,...,a,adiagD2211U,LAbUxxLD22fxBxbLDf,ULDB12122B221fxBxkk称为高斯—塞德尔迭代法(公式),用量表示的形式为(3)SOR迭代(4)交替方向迭代法（PR法）迭代格式为：对于单参数PR法，对于多参数，(5)共轭梯度法算法步骤如下：［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］,...,,k,n,,ixaxabaxijnij)k(jij)k(jijiii)k(i21021111111))1(()(1DRLDTb)(1LDdhucos)11/(22opt2121,,1,1,1,,122LLLLuuuLuuujijijijijijibuLIuLIbuLIuLIkkkkkkkkkk211122211)()()()(hhoptsin222sina....2,1)11(421k221hkahk其中（１）如果，终止算法，输出；否则下行；（２）计算：（３）计算：（４）置，转入（１）．(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;whileflagu_1=B*u_1_0+g;ifnorm(u_1-u_1_0,inf)errflag=0;endu_1_0=u_1;m_1=m_1+1;%迭代次数值enduu=zeros(n+1);formm=1:n-1fornn=1:n-1uu(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迭代法function[u_2,m_2]=Gauss_Seidel_Solve(A,b,n,err)ifnargin==2err=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;whileflagu_2=R*u_2_0+g;ifnorm(u_2-u_2_0,inf)errflag=0;endu_2_0=u_2;m_2=m_2+1;%迭代次数值enduu=zeros(n+1);formm=1:n-1fornn=1:n-1uu(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迭代法function[u_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;whileflagu_3=B1*u_3_0+g;ifnorm(u_3-u_3_0,inf)errflag=0;endu_3_0=u_3;m_3=m_3+1;%迭代次数值enduu=zeros(n+1);formm=1:n-1fornn=1:n-1uu(nn+1,mm+1)=u_3((mm-1)*(n-1)+nn,1);endend%SOR迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];figure(4)mesh(x,y,uu);title('SOR迭代差分解图像')PR方法clc;N=16;pi=3.1416;tao_opt=0.5*(1/sin(pi/N))%L1矩阵的定义L_1=zeros((N-1)*(N-1));fori=1:(N-1)*(N-1)-1L_1(i,i)=2;ifmod(i,N-1)~=0L_1(i,i+1)=-1;L_1(i+1,i)=-1;endendL_1((N-1)^2,(N-1)^2)=2;%L2矩阵的定义L_2=zeros((N-1)*(N-1));fori=1:(N-1)*(N-1)-1L_2(i,i)=2;endL_2((N-1)^2,(N-1)^2)=2;fori=N:(N-1)^2L_2(i,i-(N-1))=-1;L_2(i-(N-1),i)=-1;end%右端项的定义b=zeros((N-1)^2,1);forj=1:N-1fori=1:N-1b((j-1)*(N-1)+i)=-(2*pi^2*exp(pi*(i/N+j/N))*sin(pi*(i/N+j/N)))/(N^2);endendI=eye((N-1)*(N-1));%PR格式的建立u=zeros((N-1)^2,1);fori=1:(N-1)^2u(i)=1;endu1=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);count=1;whilenorm(u2-u,inf)=0.01count=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);forj=1:N-1fori=1:N-1z((j-1)*(N-1)+i)=exp(pi*(i/N+j/N))*(sin(pi*i/N)*sin(pi*j/N));endendzcountuu=zeros(N+1);formm=1:N-1fornn=1:N-1uu(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))')%画图数值解共轭梯度法function[u_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;whileflaga0=((r0)'*r0)/((p0)'*(A*p0));u_5=u_5_0+a0*p0;ifnorm(u_5-u_5_0,inf)errflag=0;endu_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);formm=1:n-1fornn=1:n-1uu(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('共轭梯度法迭代差分解图像')预处理共轭梯度法function[u_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+D\D;%B=L*D*(L)';B=0.7*eye((n-1)^2)p0=b-A*u_6_0;r0=p0;s0=B\r0;flag=1;whileflaga0=((r0)'*s0)/((p0)'*(A*p0));u_6=u_6_0+a0*p0;ifnorm(u_6-u_6_0,inf)errflag=0;endu_6_0=u_6;r1=r0-a0*A*p0;s1=B\r1;b0=((r1)'*s1)/((r0)'*s0);p0=s1+b0*p0;r0