Matlab-Self-adjoint elliptic PDE代码解析

pdc31czy

已于 2024-02-26 16:43:32 修改

阅读量176

点赞数 1

分类专栏：个人笔记 MATLAB 文章标签： matlab 开发语言算法

于 2022-12-19 17:27:07 首次发布

本文链接：https://blog.csdn.net/pdc31czy/article/details/128375126

版权

个人笔记同时被 2 个专栏收录

19 篇文章 3 订阅

订阅专栏

MATLAB

7 篇文章 1 订阅

订阅专栏

1. PDE

$\nabla \cdot(p(x, y) \nabla u)-q(x, y) u=f(x, y), \quad x \in(-\pi / 2, \pi / 2), \quad y \in(\pi / 2, \pi), \\ p(x, y)=x+y \\ q(x, y)=x^2+y^2 \\ f(x, y)=\cos x \cos y-\sin x \sin y-\cos x \sin y\left(x^2+y^2\right)-2 \cos x \sin y(x+y), \\ \frac{\partial u}{\partial x}=\sin y, \quad \text { on } \quad x=-\pi / 2, \\ u=0, \quad \text { on } \quad x=\pi / 2 \\ u=\cos x, \quad \text { on } \quad y=\pi / 2 \\ u=0, \quad \text { on } \quad y=\pi$
Remark: exact solution is $\cos x \sin y$

2. FDM

$\frac{p_{i+\frac{1}{2}, j} U_{i+1, j}-\left(p_{i+\frac{1}{2}, j}+p_{i-\frac{1}{2}, j}\right) U_{i j}+p_{i-\frac{1}{2}, j} U_{i-1, j}}{h_x^2} \\ +\frac{p_{i, j+\frac{1}{2}} U_{i, j+1}-\left(p_{i, j+\frac{1}{2}}+p_{i, j-\frac{1}{2}}\right) U_{i j}+p_{i, j-\frac{1}{2}} U_{i, j-1}}{h_y^2} -q_{i j} U_{i j}=f_{i j}$
也可以写成
$\frac{p_{i+\frac{1}{2}, j}}{h_x^2} U_{i+1, j} + \frac{p_{i-\frac{1}{2}, j}}{h_x^2} U_{i-1, j} + \frac{p_{i, j+\frac{1}{2}}}{h_y^2} U_{i, j+1} + \frac{p_{i, j-\frac{1}{2}}}{h_y^2} U_{i, j-1} \\+ \left(- \frac{\left(p_{i+\frac{1}{2}, j}+p_{i-\frac{1}{2}, j}\right)}{h_x^2} -\frac{\left(p_{i, j+\frac{1}{2}}+p_{i, j-\frac{1}{2}}\right)}{h_y^2} -q_{i j}\right) U_{i j}=f_{i j}$

3.Matlab Code

初始化

clearvars;
close all;
clc;

a = -0.5*pi; b= 0.5*pi; c = pi/2; d=pi;

m=128; n=128;

hx = (b-a)/m; hx1 = hx*hx; x=zeros(m+1,1);
for i=1:m+1
  x(i) = a + (i-1)*hx;
end

hy = (d-c)/n; hy1 = hy*hy; y=zeros(n+1,1);
for i=1:n+1
  y(i) = c + (i-1)*hy;
end

M = (n-1)*m; A = sparse(M,M); bf = zeros(M,1);

在这里插入图片描述

for j = 1:n-1
  for i=1:m
    k = i + (j-1)*m;
    bf(k) = f(x(i),y(j+1));
    p1=p(x(i)+hx/2,y(j+1))+p(x(i)-hx/2,y(j+1));
    p2=p(x(i),y(j+1)+hy/2)+p(x(i),y(j+1)-hy/2);
    A(k,k) = -p1/hx1-p2/hy1-q(x(i),y(j+1));

对于left boundary condition（Neumann）：

在这里插入图片描述

if i == 1
        A(k,k+1) = p1/hx1;
        bf(k) = bf(k) + 2*p(x(i)-hx/2,y(j+1))*ux(y(j+1))/hx;

对于right boundary condition：

在这里插入图片描述

 else
       if i==m
         A(k,k-1) = p(x(i)-hx/2,y(j+1))/hx1;
         bf(k) = bf(k) - p(x(i)+hx/2,y(j+1))*ue(x(i+1),y(j+1))/hx1;

对于x方向（lb和rb除外）：

在这里插入图片描述

       else
	  A(k,k-1) = p(x(i)-hx/2,y(j+1))/hx1;
	  A(k,k+1) = p(x(i)+hx/2,y(j+1))/hx1;
       end
    end

对于lower boundary condition:

在这里插入图片描述

if j == 1
        A(k,k+m) = p(x(i),y(j+1)+hy/2)/hy1;
         bf(k) = bf(k) -p(x(i),y(j+1)-hy/2)*ue(x(i),c)/hy1;

对于upper boundary condition:

在这里插入图片描述

else
       if j==n-1
         A(k,k-m) = p(x(i),y(j+1)-hy/2)/hy1;
          bf(k) = bf(k) -p(x(i),y(j+1)+hy/2)*ue(x(i),d)/hy1;

对于y方向（lowerb和ub除外）：

在这里插入图片描述

else
          A(k,k-m) = p(x(i),y(j+1)-hy/2)/hy1;
          A(k,k+m) = p(x(i),y(j+1)+hy/2)/hy1;
       end
     end

  end
end

利用AU = F 求解

U = A \bf;

画图

%--- Transform back to (i,j) form to plot the solution ---


j = 1;
for k=1:M
  i = k - (j-1)*m ;
  u(i,j) = U(k);
  u2(i,j) = ue(x(i),y(j+1));
  j = fix(k/m) + 1;
end

% Analyze abd Visualize the result.

e = max( max( abs(u-u2)))        % The maximum error
x1=x(1:m); y1=y(2:n);

mesh(y1,x1,u);
title('The solution plot');
xlabel('y');
ylabel('x');
figure(2);
mesh(y1,x1,u-u2);
title('The error plot');
xlabel('y');
ylabel('x');