2.2 Neumann境界値問題

poisson1n_nh.m

% poisson1n_nh.m
% 1次元領域におけるPoisson -△u=f (非同次Neumann境界条件) を差分法で解く
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson1n_nh.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/poisson2n-nonhomo.pdf
% (2024/8/27)。
function poisson1n_nm(N)
  arguments
    N=4 % デバッグ用
  end
  function K = poisson1n_mat(W,N)
    h=W/N;
    I=speye(N+1,N+1);
    v=ones(N,1);
    J=sparse(diag(v,1)+diag(v,-1));
    K=2*I-J;
    K(1,1)=1; K(1,2)=0; K(2,1)=0;
    K(N+1,N+1)=1; K(N+1,N+1)=1;
    K=1/h^2*K;
  end
  % 想定している厳密解
  function val = exact_solution(x)% =x*(1-x)=x-x^2
    val = x .* (1 - x);
  end
  % f=-△u=-d^2u/dx^2
  function val = f(x)% = -u''(x) = 2
    val = 0 .* x + 2.0;
  end
  % ∂u/∂n(a)=-∂u/∂x(a)
  function val = b_l(x)% =-u'(0)=-(1-2x)|x=0
    val = -1;
  end
  % ∂u/∂n(b)=∂u/∂x(b)
  function val = b_r(x)% 1-2x|x=1
    val = -1;
  end
  % 領域
  a=0; b=1;
  beta=-1;
  h=(b-a)/N;
  fprintf('N=%d, h=%g\n', N, h);
  A=poisson1n_mat(b-a,N);
  if N <= 10
    disp('A='); disp(full(A));
  end
  fprintf('det(A)=%e\n', det(A))
  fprintf('cond(A)=%e\n', condest(A));

  X=linspace(a,b,N+1);
  F=f(X');
  if N <= 10
    disp('X='); disp(X);
    disp('F=f(X)='); disp(F');
  end
  fprintf('sum=%g\n', sum(F(:,1)));
  fprintf('beta=%g, f(a)/2=%g, f(b)/2=%g, h=%g\n',...
          b_r(b), f(a)/2, f(b)/2, h);
  integral=h*(f(a)/2+sum(F(2:N,1))+f(b)/2);
  alpha=-(integral+b_r(b));
  fprintf('積分=%g, 整合条件の成り立つα=-u''(0)=-(%g+(%g))=%g, 想定値=%g\n',...
          integral, integral, b_r(b), alpha, b_l(a));
  F(1)=0.0; F(N+1)=0.5*F(N+1)+b_r/h;
  if N <= 10
    fprintf('F= (F(%d)=f(1)/2+beta/h=%g/2+(%g/%g)=%g)\n', N+1, f(b), beta, h, F(N+1));
    disp(F);
  end
  %whos A F
  u=A\F;
  plot(X,u); title('差分解(左端で0のもの)のグラフ'); drawnow; shg;
  error=norm(u-exact_solution(X'),inf);
  fprintf('max norm of error=%e\n', error);
end

入手するにはターミナルで

curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson1n_nh.m

これを実行するには、例えば

cp -p poisson1n_nh.m ~/Documents/MATLAB

とコピーして (私はシンボリック・リンクをはることにしている)、 MATLAB の中で

>> poisson1n_nh

(分割数はデフォールトの $N=4$ が採用される。)

あるいは、分割数 $N$ を指定するため

>> poisson1n_nh(10)

のようにする。