桂田 [4] という解説を書いた。
poisson2n_nh.m |
% poisson2n_nh.m % 長方形領域におけるPoisson -△u=f (非同次Neumann境界条件) を差分法で解く % 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2n_nh.m % 一応の完成をして、poisson2n_nh.m として公開した (2024/8/27)。 function poisson2n_nh(Nx,Ny) arguments Nx=50 Ny=50 end % Neumann境界条件の場合の係数行列 % poisson2n_mat.m function A = poisson2n_mat(W,H,Nx,Ny) Ix=speye(Nx+1,Nx+1); Iy=speye(Ny+1,Ny+1); vx=ones(Nx,1); Jx=sparse(diag(vx,1)+diag(vx,-1)); Jx(1,2)=2; Jx(Nx+1,Nx)=2; vy=ones(Ny,1); Jy=sparse(diag(vy,1)+diag(vy,-1)); Jy(1,2)=2; Jy(Ny+1,Ny)=2; hx = W / Nx; hy = H / Ny; betax = 1.0 / hx^2; betay = 1.0 / hy^2; Kx=betax * (2*Ix - Jx); Ky=betay * (2*Iy - Jy); % row major (y changes first, l=j+(Ny+1)*i) A=kron(Kx,Iy) + kron(Ix,Ky); % column major (x chandes first, l=i+(Nx+1)*j) % A=kron(Iy,Kx) + kron(Ky,Ix); end % 長方形領域 (a,b)×(c,d) a=0; b=1; c=0; d=1; % function u = exact_solution(x, y) u = sin(pi * x) .* sin(pi * y); end % Poisson 方程式の右辺 (exact_solution の -△) function val = f(x, y) val = 2 * pi^2 * sin(pi * x) .* sin(pi * y); end % 長方形の下の辺での b() = -∂u/∂y=-π sin πx cos π0=-π sin πx function val = b_b(x) val = - pi * sin(pi * x); end % 長方形の上の辺での b() = ∂u/∂y=π sin πx cos π・1=-π sin πx function val = b_t(x) val = - pi * sin(pi * x); end % 長方形の左の辺での b() = -∂u/∂x=-π cos π0 sin πy=-π sin πy function val = b_l(y) val = - pi * sin(pi * y); end % 長方形の右の辺での b() = ∂u/∂x=π cosπ・1 sin πy=-π sin πy function val = b_r(y) val = - pi * sin(pi * y); end % 誤差を測るためのノルム (最大値ノルム) function val = supnorm(a) [m,n]=size(a); val = norm(reshape(a,m*n,1),Inf); end % 格子への分割 %Nx=100; %Ny=100; hx=(b-a)/Nx; hy=(d-c)/Ny; betax=1.0/hx^2; betay=1.0/hy^2; disp(sprintf('Nx=%d, Ny=%d, hx=%f, hy=%f, betax=%f, betay=%f', ... Nx, Ny, hx, hy, betax, betay)); % 差分方程式 A U^{n+1}=B U^n の行列 A=poisson2n_mat(b-a,d-c,Nx,Ny); % 対称にするためのスケーリング用の行列 Sx=sparse(diag([1/2;ones(Nx-1,1);1/2])); Sy=sparse(diag([1/2;ones(Ny-1,1);1/2])); S=kron(Sy,Sx); % Aを対称行列に直す A=S*A; % 対称かどうかチェックする er=A-A'; if norm(er,inf)>=1 disp('Aが対称でない。') pause else %disp('Aは対称である。'); end clear er; N=(Nx+1)*(Ny+1); if N <= 50 disp('A='); disp(full(A)); end % 格子点の座標ベクトル x=(x_1,x_2,...,x_{Nx+1}), y=(y_1,y_2,...,y_{Ny+1}) X=linspace(a,b,Nx+1)'; Y=linspace(c,d,Ny+1)'; % 格子点のx,y座標の配列 X={X_{ij}}, Y={Y_{ij}} [x,y]=meshgrid(X,Y); % Poisson方程式の右辺の関数値 (連立1次方程式に必要) u=f(x,y); close all hidden meshc(x,y,u); axis([a b c d -50 50]); title('graph of f'); drawnow; shg % 計算に用いるため、f をベクトル U とする。 U=reshape(u,N,1); if N <= 50 disp('f='); disp(u); end % 境界値 (左の辺、右の辺) blvector=b_l(Y); brvector=b_r(Y); lindex=1:Ny+1; rindex=Nx*(Ny+1)+1:(Nx+1)*(Ny+1); if Ny+1 <= 20 disp('blvector='); disp(blvector); disp('brvector='); disp(brvector); disp('lindex='); disp(lindex); disp('rindex='); disp(rindex); end % 境界値 (下の辺、上の辺) bbvector=b_b(X); btvector=b_t(X); bindex=1:Ny+1:(Nx+1)*(Ny+1); tindex=Ny+1:Ny+1:(Nx+1)*(Ny+1); if Nx+1 <= 20 disp('bbvector='); disp(bbvector); disp('btvector='); disp(btvector); disp('bindex='); disp(bindex); disp('tindex='); disp(tindex); end % 境界条件により分かっている値を右辺に移項 U(lindex)=U(lindex)+2*hx*betax*blvector; U(rindex)=U(rindex)+2*hx*betax*brvector; U(bindex)=U(bindex)+2*hy*betay*bbvector; U(tindex)=U(tindex)+2*hy*betay*btvector; % 係数行列を対称化するためのスケーリングを U にも施す U=S*U; if N <= 50 disp('右辺のベクトル='); disp(reshape(U,Ny+1,Nx+1)); end % U_{0,0}=0 として解く A(1,1)=1; A(1,2)=0; A(1,Ny+2)=0; A(2,1)=0; A(Ny+2,1)=0; U(1)=0.0; if N <= 20 disp('A='); disp(full(A)); end % LU分解 [AL,AU,ap]=lu(A,'vector'); if N<= 20 disp('U='); disp(U); end % 連立1次方程式を解いて差分解を求める U=AU\(AL\U(ap,:)); % 差分解を meshgrid に収める u=reshape(U,Ny+1,Nx+1); % if N <= 50 disp('u='); disp(u); end diff=u(1,1)-exact_solution(1,1); %fprintf('差分解(0,0)-厳密解(0,0)=%g\n', diff); u=u-diff; % 解のグラフ minu=min(U); maxu=max(U); figure; meshc(x,y,u); axis([a b c d -1 1]); title('graph of the approximate solution'); drawnow; shg; % 誤差を調べる error = supnorm(u-exact_solution(x,y)); disp(sprintf('Nx=%d, Ny=%d, ||error||=%e, ||u||=%e',... Nx, Ny, error, supnorm(u))); % 比較のため厳密解を表示する meshc(x,y,exact_solution(x,y)); axis([a b c d -1 1]); title('graph of the exact solution'); drawnow; shg end |
入手するには、ターミナルで
curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2n_nh.m |
cp -p poisson2n_nh.m ~/Documents/MATLAB |
>> poisson2n_nh(分割数はデフォールトの , が採用される。) |
>> poisson2n_nh(100,100) |