next up previous contents
Next: 2.2 Neumann境界値問題 Up: 2 1次元Poisson方程式 Previous: 2 1次元Poisson方程式

2.1 Dirichlet境界値問題

$\displaystyle -u''(x)=f(x)$   ( $ x\in (0,L)$)$\displaystyle ,\quad u(0)=u(L)=0 $

$\displaystyle h=\frac{L}{N},\quad x_i=ih$   ( $ i=0,1,\cdots,N$)$\displaystyle . $

$\displaystyle \frac{1}{h^2} \begin{pmatrix} 2 & -1 \\ -1 & 2 & -1 \\ & -1... ...begin{pmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_{N-1}) \end{pmatrix}. $

laplacian1d.m

% laplacian1d.m
% curl -O https://m-katsurada.sakura.ne.jp/program/fdm/laplacian1d.m
% -(d/dx)^2 の差分近似
% [0,L] を N 等分した時の行列は h=L/N; a=laplacian1d(N-1)/(h*h);
function a=laplacian1d(m)
  e=ones(m,1);
  a=spdiags([-e 2*e -e],-1:1,m,m); % 疎行列なので圧縮した形式
end


>> a=laplacian1d(3));
>> full(a)

ans =
     2    -1     0
    -1     2    -1
     0    -1     2

$ I=(0,L)$ を $ N$ 等分する差分近似の場合は 
>> L=1;
>> n=5;
>> h=L/n;
>> a=laplacian1d(n-1)/(h*h);
>> full(a)

ans =

   50.0000  -25.0000         0         0
  -25.0000   50.0000  -25.0000         0
         0  -25.0000   50.0000  -25.0000
         0         0  -25.0000   50.0000

$ I=(0,L)$ における

$\displaystyle -u''(x)=1,\quad u(0)=u(L)=0 $

を解くには ($ f\equiv 1$ とするのは安直だけど)
poisson1d_test1.m

L=1
n=10
h=L/n
x=linspace(0,L,n+1);
a=laplacian1d(n-1)/(h*h);
f=ones(n-1,1);
u=a\f;
u=[0; u; 0];
plot(x,u)
figure(gcf)
とする。

$ f(x)=\sin x$ の場合は
poisson1d_test2.m

L=1
n=10
h=L/n
x=linspace(0,L,n+1);
a=laplacian1d(n-1)/(h*h);
f=sin(x(2:n))';
u=a\f;
u=[0; u; 0];
plot(x,u)
figure(gcf)
とする。

$ f(x)=x(1-x)$ の場合は
poisson1d_test3.m

L=1
n=10
h=L/n
x=linspace(0,L,n+1);
a=laplacian1d(n-1)/(h*h);
f=x.*(1-x);
f=f(2:n)';
u=a\f;
u=[0; u; 0];
plot(x,u)
figure(gcf)
とする。


(2024/9 改訂版) 非同次Dirichlet境界値問題のプログラム。

poisson1d_nh.m

% poisson1d_nh.m
% 1次元領域におけるPoisson -△u=f (非同次Dirichlet境界条件) を差分法で解く
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson1d_nh.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/poisson1d.pdf
% (2024/9/4)。
function poisson1d_nh(N)
  arguments
    N=4 % デバッグ用
  end
  function K = poisson1d_mat(W,N)
    h=W/N;
    I=speye(N-1,N-1);
    v=ones(N-2,1);
    J=sparse(diag(v,1)+diag(v,-1));
    K=2*I-J;
    K=1/h^2*K;
  end
  % 想定している厳密解
  function val = exact_solution(x)% =x^4
    val = x .^ 4;
    %val = (x - 1.0/2).^2;
  end
  % f=-△u=-d^2u/dx^2
  function val = f(x)% = -u''(x) = -12x^2
    val = - 12 * x.^2;
    %val = 0 .* x - 2.0;
  end
  % 領域
  a=0; b=1;
  % 境界値
  alpha = exact_solution(a);
  beta = exact_solution(b);
  fprintf('a=%f, b=%f, u(a)=%f, u(b)=%f\n', a, b, alpha, beta);
  % 差分近似の行列
  h=(b-a)/N;
  fprintf('N=%d, h=%g\n', N, h);
  A=poisson1d_mat(b-a,N);
  if N <= 10
    disp('A='); disp(full(A));
  end
  % 連立方程式の右辺
  X=linspace(a,b,N+1);
  F=f(X(2:N)');
  F(1)=F(1)+alpha/h^2;
  F(N-1)=F(N-1)+beta/h^2;
  if N <= 10
    disp('X='); disp(X);
    disp('F=f(X)='); disp(F');
  end
  % 差分解を求める
  u=zeros(N+1,1);
  u(2:N)=A\F;
  % 境界値
  u(1)=alpha; u(N+1)=beta;
  % グラフを描く
  plot(X,u); title('差分解のグラフ'); drawnow; shg;
  figure
  plot(X,exact_solution(X')); title('厳密解のグラフ'); drawnow; shg;
  exact=exact_solution(X');
  error=norm(u-exact,inf);
  fprintf('max norm of error=%e\n', error);
end

入手するにはターミナルで 
curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson1d_nh.m

これを実行するには、例えば 
cp -p poisson1d_nh.m ~/Documents/MATLAB
とコピーして (私はシンボリック・リンクをはることにしている)、 MATLAB の中で 
>> poisson1d_nh
(分割数はデフォールトの $ N=4$ が採用される。)
あるいは、分割数 $ N$ を指定するため 
>> poisson1d_nh(10)
のようにする。

上のサンプル・プログラムでは、 $ u(x)=x^4$ と選び、 $ f(x)=-u''(x)=-12x^2$, $ \alpha=u(0)=0$, $ \beta=u(1)=1$ としてある。 $ N$ を $ 10$ から倍々にしていくと以下のようになった。 
N=10, h=0.1
max norm of error=2.500000e-03

N=20, h=0.05
max norm of error=6.250000e-04

N=40, h=0.025
max norm of error=1.562500e-04

N=80, h=0.0125
max norm of error=3.906250e-05

N=160, h=0.00625
max norm of error=9.765625e-06

N=320, h=0.003125
max norm of error=2.441406e-06
$ N$が$ 2$倍になると誤差は $ \frac{1}{4}$ 倍になる。 誤差は $ O(1/N^2)$ になっているようで、ひとまずは納得できる。


桂田 祐史