请选择 进入手机版 | 继续访问电脑版

Poisson方程的五点差分格式例题求解-Matlab实现

[复制链接]
太阳神鹰 发表于 2021-1-3 12:15:31 | 显示全部楼层 |阅读模式 打印 上一主题 下一主题
文章目次



1.前言

本文使用Matalab软件,将应用偏微分方程数值解中椭圆形方程的五点差分格式求解一道Poisson方程的第一边值例题,通过五点差分格式及边值条件得到相应的差分方程组                              K                      U                      =                      F                          KU=F               KU=F后,接纳Gauss-seidel迭代法对其求解即可得到数值解,并将数值解与真解作比力.
此中五点差分格式将直接给出,其构造过程从略.
2.例题

构造Poisson方程第一边值如下:

接纳五点差分格式求解并与剖析解                               u                      (                      x                      ,                      y                      )                      =                               x                         2                                       y                         2                                  u(x,y)=x^2y^2               u(x,y)=x2y2举行比力.
(                              x                      和                      y                          x和y               x和y方向均                              n                          n               n平分,取                                       h                         1                              =                               h                         2                              =                      h                      =                      1                      /                      n                          h_1=h_2=h=1/n               h1​=h2​=h=1/n,                              x                      和                      y                          x和y               x和y方向上的节点序号分别用                              i                      ,                      j                          i,j               i,j体现)
3.符号说明

表1:
符号说明                                             n                                      n                        n区间平分数                                             h                                      h                        h区间剖分步长                                             K                                      K                        K方程组系数矩阵                                             U                                      U                        U方程组未知向量                                             F                                      F                        F方程组右端项                                             u                                      u                        u方程组数值解向量                                                         U                                  1                                                 U_1                        U1​将                                             u                                      u                        u的元素放到矩阵                                                         U                                  1                                                 U_1                        U1​中                                                         U                                  2                                                 U_2                        U2​剖析解矩阵                                             N                                      N                        N迭代次数上限                                             e                               p                                      ep                        ep迭代误差限                                             k                                      k                        k迭代次数                                             Q                                      Q                        Q区域界限节点值的算术均匀值4.思路

主要思路是将二维问题当成一维问题求解,即二维的                              (                      n                      +                      1                               )                         2                                  (n+1)^2               (n+1)2个节点按顺序分行放到                              (                      n                      +                      1                               )                         2                                  (n+1)^2               (n+1)2维向量                              U                          U               U中.
主要步调如下:
1.确定步长后对区域举行网格匀称剖分(                              x                      ,                      y                          x,y               x,y方向的步长                              h                          h               h相等);
2.所求方程组为                              K                      U                      =                      F                          KU=F               KU=F,根据五点差分格式分别给                              K                          K               K和                              F                          F               F赋值(矩阵中包罗界限点);
(说明:                              U                          U               U是一个                              (                      n                      +                      1                               )                         2                                  (n+1)^2               (n+1)2维向量,即将求解区域中所有节点放到一个向量                              U                          U               U中,当作一维问题举行求解;其次,将界限点放入方程组中可以使                              K                          K               K具有更好的性质,更好地防止求解中不收敛的情况,处置惩罚界限点时将其当成未知的即可)
3.由界限点值的算术均匀值作为                              U                          U               U的初始值,以淘汰迭代次数;
4.给定迭代次数上限                              N                          N               N和误差限,由                              G                      a                      u                      s                      s                      −                      s                      e                      i                      d                      e                      l                          Gauss-seidel               Gauss−seidel迭代求解方程组                              K                      U                      =                      F                          KU=F               KU=F得到数值解向量                              u                          u               u;
5.为方便观察,将解向量                              u                          u               u的元素放入                              (                      n                      +                      1                      )                      ∗                      (                      n                      +                      1                      )                          (n+1)*(n+1)               (n+1)∗(n+1)阶矩阵                                       U                         1                                  U_1               U1​中,并与由剖析解                               u                      (                      x                      ,                      y                      )                      =                               x                         2                                       y                         2                                  u(x,y)=x^2y^2               u(x,y)=x2y2产生的解矩阵                                       U                         2                                  U_2               U2​比力;
5.求解步调

分别网格数                              n                          n               n选取为7.
1.右端项为                                       f                                   i                            j                                       =                      −                      2                      [                      (                      i                      h                               )                         2                              +                      (                      j                      h                               )                         2                              ]                          f_{ij}=-2[(ih)^2+(jh)^2]               fij​=−2[(ih)2+(jh)2],则得到方程的五点差分格式为

相应的                              G                      a                      u                      s                      s                      −                      s                      e                      i                      d                      e                      l                          Gauss-seidel               Gauss−seidel迭代公式为

2.构造求解矩阵时,对特殊范例的节点应根据界限条件将格式适当化简:
(1)左界限点                                       u                                   0                            j                                       (                      j                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                      4                               u                                   0                            j                                       =                      0                          u_{0j}(j=0,1,...,n):4u_{0j}=0               u0j​(j=0,1,...,n):4u0j​=0
(2)下界限点                                       u                                   i                            0                                       (                      i                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                      4                               u                                   i                            0                                       =                      0                          u_{i0}(i=0,1,...,n):4u_{i0}=0               ui0​(i=0,1,...,n):4ui0​=0
(3)右界限点                                       u                                   n                            j                                       (                      j                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                      4                               u                                   n                            j                                       =                      4                      (                      j                      h                               )                         2                                  u_{nj}(j=0,1,...,n):4u_{nj}=4(jh)^2               unj​(j=0,1,...,n):4unj​=4(jh)2
(4)上界限点                                       u                                   i                            n                                       (                      i                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                      4                               u                                   i                            n                                       =                      4                      (                      i                      h                               )                         2                                  u_{in}(i=0,1,...,n):4u_{in}=4(ih)^2               uin​(i=0,1,...,n):4uin​=4(ih)2
(5)左界限点的相邻内点                                       u                                   1                            j                                       (                      j                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                          u_{1j}(j=0,1,...,n):               u1j​(j=0,1,...,n):
                               −                      (                               u                                   2                            j                                       −                      2                               u                                   1                            j                                       )                      −                      (                               u                                   1                            ,                            j                            +                            1                                       −                      2                               u                                   1                            j                                       +                               u                                   1                            ,                            j                            −                            1                                       )                      =                               h                         2                                       f                                   1                            j                                           -(u_{2j}-2u_{1j})-(u_{1,j+1}-2u_{1j}+u_{1,j-1})=h^2f_{1j}               −(u2j​−2u1j​)−(u1,j+1​−2u1j​+u1,j−1​)=h2f1j​
(6)下界限点的相邻内点                                       u                                   i                            1                                       (                      i                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                          u_{i1}(i=0,1,...,n):               ui1​(i=0,1,...,n):
                               −                      (                               u                                   i                            +                            1                            ,                            1                                       −                      2                               u                                   i                            1                                       +                               u                                   i                            −                            1                            ,                            1                                       )                      −                      (                               u                                   i                            2                                       −                      2                               u                                   i                            1                                       )                      =                               h                         2                                       f                                   i                            1                                           -(u_{i+1,1}-2u_{i1}+u_{i-1,1})-(u_{i2}-2u_{i1})=h^2f_{i1}               −(ui+1,1​−2ui1​+ui−1,1​)−(ui2​−2ui1​)=h2fi1​
(7)右界限点的相邻内点                                       u                                   n                            −                            1                            ,                            j                                       (                      j                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                          u_{n-1,j}(j=0,1,...,n):               un−1,j​(j=0,1,...,n):
                               −                      (                      −                      2                               u                                   n                            −                            1                            ,                            j                                       +                               u                                   n                            −                            2                            ,                            j                                       )                      −                      (                               u                                   n                            −                            1                            ,                            j                            +                            1                                       −                      2                               u                                   n                            −                            1                            ,                            j                                       +                               u                                   n                            −                            1                            ,                            j                            −                            1                                       )                      =                               h                         2                                       f                                   n                            −                            1                            ,                            j                                       +                      (                      j                      h                               )                         2                                  -(-2u_{n-1,j}+u_{n-2,j})-(u_{n-1,j+1}-2u_{n-1,j}+u_{n-1,j-1})=h^2f_{n-1,j}+(jh)^2               −(−2un−1,j​+un−2,j​)−(un−1,j+1​−2un−1,j​+un−1,j−1​)=h2fn−1,j​+(jh)2
(8)上界限点的相邻内点                                       u                                   i                            ,                            n                            −                            1                                       (                      i                      =                      0                      ,                      1                      ,                      .                      .                      .                      ,                      n                      )                      :                          u_{i,n-1}(i=0,1,...,n):               ui,n−1​(i=0,1,...,n):
                               −                      (                               u                                   i                            +                            1                            ,                            n                            −                            1                                       −                      2                               u                                   i                            ,                            n                            −                            1                                       +                               u                                   i                            −                            1                            ,                            n                            −                            1                                       )                      −                      (                      −                      2                               u                                   i                            ,                            n                            −                            1                                       +                               u                                   i                            ,                            n                            −                            2                                       )                      =                               h                         2                                       f                                   i                            ,                            n                            −                            1                                       +                      (                      i                      h                               )                         2                                  -(u_{i+1,n-1}-2u_{i,n-1}+u_{i-1,n-1})-(-2u_{i,n-1}+u_{i,n-2})=h^2f_{i,n-1}+(ih)^2               −(ui+1,n−1​−2ui,n−1​+ui−1,n−1​)−(−2ui,n−1​+ui,n−2​)=h2fi,n−1​+(ih)2
对于同时满足上述多个情况的节点,只需根据实际对上述相应格式举行一定的叠加即可.
3.根据五点差分格式和界限条件(步调2)对                              K                      ,                      U                      和                      F                          K,U和F               K,U和F赋值(此中选取界限值的算术均匀值                              Q                      =                      0.3242                          Q=0.3242               Q=0.3242为初始近似值赋于未知向量                              U                          U               U.)
4.取最大迭代次数N=500,迭代误差限                               e                      p                          ep               ep分别取                              1                      e                      −                      2                          1e-2               1e−2和                              1                      e                      −                      4                          1e-4               1e−4,接纳                              G                      a                      u                      s                      s                      −                      s                      e                      i                      d                      e                      l                          Gauss-seidel               Gauss−seidel迭代求解方程组                              K                      U                      =                      F                          KU=F               KU=F,将求解效果放到                              (                      n                      +                      1                               )                         2                                  (n+1)^2               (n+1)2维向量                              u                          u               u中.
5.将                              u                          u               u中的元素按一定序次存储到                              (                      n                      +                      1                      )                          (n+1)               (n+1)阶矩阵                                       U                         1                                  U_1               U1​中,同时通过理论解                              u                      (                      x                      ,                      y                      )                      =                               x                         2                                       y                         2                                  u(x,y)=x^2y^2               u(x,y)=x2y2产生相应节点的函数值并存储到                              (                      n                      +                      1                      )                          (n+1)               (n+1)阶矩阵                                       U                         2                                  U_2               U2​中,将                                       U                         1                              ,                               U                         2                                  U_1,U_2               U1​,U2​的元素举行比力和分析.
6.求解效果

表2:理论解得到的节点函数值(                                       U                         2                              中                      元                      素                          U_2中元素               U2​中元素)
                                                          u                                               i                                     j                                                             u_{ij}                        uij​                                             i                               =                                      i=                        i= 01234567                                             j                               =                               0                                      j=0                        j=000000000100.00040.00170.00370.00670.01040.01500.0204200.00170.00670.01500.02670.04160.06000.0816300.00370.01500.03370.06000.09370.13490.1837400.00670.02670.06000.10660.16660.23990.3265500.01040.04160.09370.16660.26030.37480.5102600.01500.06000.13490.23990.37480.53980.7347700.02040.08160.18370.32650.51020.73471.0000表3:                              e                      p                      =                      1                      e                      −                      2                          ep=1e-2               ep=1e−2时得到的节点函数值(                                       U                         1                              中                      元                      素                          U_1中元素               U1​中元素,此时迭代次数                              k                      =                      8                          k=8               k=8)
                                                          u                                               i                                     j                                                             u_{ij}                        uij​                                             i                               =                                      i=                        i= 01234567                                             j                               =                               0                                      j=0                        j=000000000100.01250.02030.02350.02330.02190.02050.0204200.02030.03540.04540.05250.05960.06870.0816300.02350.04540.06600.08770.11310.14450.1837400.02330.05250.08770.13060.18350.24830.3265500.02190.05960.11310.18350.27230.38080.5102600.02050.06870.14450.24830.38080.54280.7347700.02040.08160.18370.32650.51020.73471.0000表4:                              e                      p                      =                      1                      e                      −                      4                          ep=1e-4               ep=1e−4时得到的节点函数值(                                       U                         1                              中                      元                      素                          U_1中元素               U1​中元素,此时迭代次数                              k                      =                      29                          k=29               k=29)
                                                          u                                               i                                     j                                                             u_{ij}                        uij​                                             i                               =                                      i=                        i= 01234567                                             j                               =                               0                                      j=0                        j=000000000100.00050.00190.00400.00690.01050.01510.0204200.00190.00700.01530.02700.04190.06010.0816300.00400.01530.03410.06030.09400.13510.1837400.00690.02700.06030.10690.16680.24000.3265500.01050.04190.09400.16680.26050.37490.5102600.01510.06010.13510.24000.37490.53980.7347700.02040.08160.18370.32650.51020.73471.0000通过mesh函数分别将表2-表4的元素可视化,依次得到图1,图2,图3如下.
图1:

图2:

图3:

7.总结

通过效果可看出,当迭代误差限不绝减小到时,迭代次数从8次上升到29次,方程组数值解也更靠近于理论解.本次实验选取巨细适中,可以大概更广泛地反映求解区域中各节点的数值解取值,同时也方便直观地举行比力.可以推断当误差限趋于0时,迭代次数将趋于无穷,此时五点差分格式盘算得到的数值解将无限趋近于理论解.但由于选取的误差限较小,节点个数选取仍然较小,我们不易直观地从图1至图3中直观地观察到两次数值效果同理论解之间的差别.
8.Matlab代码

1.主函数
  1. %% 所求方程为KU=F......将二维问题当成一维问题求解,即二维的(n+1)^2个节点按顺序分行放到(n+1)^2维向量U中clcn=7;       %网格数h=1/n;     %步长F=zeros((n+1)^2,1); %KU=FK=zeros(size(F));   %K为(n+1)^2阶方阵U=ones((n+1)^2,1);  %K为(n+1)^2维向量U0=zeros((n+1)^2,1);%U0用作存边值条件,为初始向量做准备for i=1:n+1    %对U0中界限点的位置赋值    U0(i)=0;            %下界限点    U0(n^2+n+i)=(i*h)^2;%右界限点    U0(1+(i-1)*(n+1))=0;%左界限点    U0(i*(n+1))=(i*h)^2;%上界限点end%% 下面临K赋值for i=1:(n+1)^2    K(i,i)=4;endfor i=1:n-1 %五点中的右点    for j=2+i*(n+1):(n-1)+i*(n+1)        K(j,j+1)=-1;    endendfor i=1:n-1 %五点中的左点    for j=3+i*(n+1):n+i*(n+1)        K(j,j-1)=-1;    endendfor i=1:n-2 %五点中的上点    for j=2+i*(n+1):n+i*(n+1)        K(j,j+n+1)=-1;    endendfor i=2:n-1 %五点中的下点    for j=2+i*(n+1):n+i*(n+1)        K(j,j-n-1)=-1;    endend%% 下面临F赋值%下面赋值非特殊的点for i=1:n-1    for j=2+i*(n+1):n+i*(n+1)        F(j)=-2*((floor(j/(n+1)))^2+(mod(j,n+1)-1)^2)*h^4;    endend%下面赋值与界限点相邻的内点(左边和下边为齐次,不消管)for i=1:n-1 %右边的点    F(n+i*(n+1))=F(n+i*(n+1))+(i*h)^2;  endfor i=1:n-1%上边的点    F((n+1)*(n-1)+i+1)=F((n+1)*(n-1)+i+1)+(i*h)^2;end%下面赋值F中界限点for i=i:n+1    F(i)=0;%下界限点endfor i=1:n-1    F(1+i*(n+1))=0;%左界限点    F((i+1)*(n+1))=4*(i*h)^2;%右界限点endfor i=(n+1)*n+1:(n+1)^2 %上界限点    if i==(n+1)^2        F(i)=4*(n^2)*(h^2);    else        F(i)=4*((mod(i,n+1)-1)^2)*(h^2);    endend%% 下面临U赋初值d=sum(sum(U0))/4/n;  %使用界限点值的算术均匀值作为U的初值以加速迭代次数for i=1:(n+1)^2    U(i)=d;end%% 使用高斯赛德迭代求解KU=Fep=1e-4;  %误差限,可根据需要更改,本次我们使用1e-2和1e-4N=500;    %最大迭代次数u=Gauss(K,F,U,ep,N);  %将求解效果存到向量u中%% 下面将解u放到矩阵U1中,并与真解U2作比力U1=zeros(n+1);%U1用于存放向量u中元素二维化后的数据U2=zeros(n+1);%U2用于存放真解for i=1:n+1%u的一维数据放到二维矩阵U1中    for j=1+(i-1)*(n+1):i*(n+1)        U1(i,j-(i-1)*(n+1))=u(j);    endendfor i=1:n+1%真解产生的节点函数值放到U2中    for j=1:n+1        U2(i,j)=(((i-1)*h)^2)*(((j-1)*h)^2);    endend%% 输出解U1U2%% 画图X=linspace(0,1,n+1);Y=linspace(0,1,n+1);mesh(X,Y,U2);%当要画数值解的图时,U2改为U1即可hold on;title('理论解图像');%当要画数值解的图时,'理论解图像'改为'数值解图像'即可
复制代码
2.Gauss-seidel迭代
[code]function x=Gauss(A,b,x0,ep,N)%用于Gauss-seidel迭代法解线性方程组Ax=b%A,b,x0分别为系数矩阵,右端向量和初始向量(初始向量默认为零向量)%ep为精度(1e-3),N为最大迭代次数(默认500次),x返回数值解向量n=length(b);if nargin
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

发布主题

专注素材教程免费分享
全国免费热线电话

18768367769

周一至周日9:00-23:00

反馈建议

27428564@qq.com 在线QQ咨询

扫描二维码关注我们

Powered by Discuz! X3.4© 2001-2013 Comsenz Inc.( 蜀ICP备2021001884号-1 )