% 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
|