非同次の問題の解説 [3] も書いた。
ここでは、 [3] で作成したプログラムを紹介する、
poisson2d_nh.m |
% poisson2d_nh.m % 長方形領域におけるPoisson -△u=f (非同次Dirichlet境界条件) を差分法で解く % 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2d_nh.m % 解説 https://m-katsurada.sakura.ne.jp/labo/text/poisson2d-nonhomo.pdf % 一応の完成をして、poisson2d_nh.m として公開した (2024/8/25)。 function poisson2d_nh(Nx,Ny) arguments Nx=50 Ny=50 end % Dirichlet境界条件の場合の係数行列 function A = poisson2d_mat(W,H,Nx,Ny) Ix=speye(Nx-1,Nx-1); Iy=speye(Ny-1,Ny-1); vx=ones(Nx-2,1); Jx=sparse(diag(vx,1)+diag(vx,-1)); vy=ones(Ny-2,1); Jy=sparse(diag(vy,1)+diag(vy,-1)); 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-1)) A=kron(Kx,Iy) + kron(Ix,Ky); % column major (x chandes first, l=i+(Nx-1)*(j-1)) % A=kron(Iy,Kx) + kron(Ky,Ix); end function A=poisson_coef(W, H, nx, ny) hx=W/nx; hy=H/ny; m=nx-1; n=ny-1; ex=ones(nx,1); ey=ones(ny,1); Lx=spdiags([-ex,2*ex,-ex],-1:1,m,m)/(hx*hx); Ly=spdiags([-ey,2*ey,-ey],-1:1,n,n)/(hy*hy); A=kron(speye(m,m),Ly)+kron(Lx,speye(n,n)); end % 長方形領域 (a,b)×(c,d) a=0; b=1; c=0; d=1; % function u = exact_solution(x, y) u = cos(pi * x) .* cos(pi * y); end % Poisson 方程式の右辺 (exact_solution の -△) function val = f(x, y) val = 2 * pi^2 * cos(pi * x) .* cos(pi * y); end % 長方形の下の辺での b() function val = b_b(x) val = exact_solution(x,c); end % 長方形の上の辺での b() function val = b_t(x) val = exact_solution(x,d); end % 長方形の左の辺での b() function val = b_l(y) val = exact_solution(a, y); end % 長方形の右の辺での b() function val = b_r(y) val = exact_solution(b, 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=poisson2d_mat(b-a,d-c,Nx,Ny); N=(Nx-1)*(Ny-1); if N <= 20 disp('A='); full(A) end %A2=poisson_coef(b-a,d-c,Nx,Ny); %err=A-A2; %disp(sprintf('A-A2=%e', supnorm(err))); %clear A2 err % LU分解 [AL,AU,ap]=lu(A,'vector'); % 格子点の座標ベクトル 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); meshc(x,y,u); axis([a b c d -50 50]); title('graph of f'); drawnow; shg % 計算に用いるため、f の内部格子点での値をベクトル U に格納する U=reshape(u(2:Ny,2:Nx),(Nx-1)*(Ny-1),1); % 境界値 (左の辺、右の辺) blvector=b_l(Y(2:Ny)); brvector=b_r(Y(2:Ny)); lindex=1:Ny-1; rindex=(Nx-2)*(Ny-1)+1:(Nx-1)*(Ny-1); % 境界値 (下の辺、上の辺) bbvector=b_b(X(2:Nx)); btvector=b_t(X(2:Nx)); bindex=1:Ny-1:(Nx-1)*(Ny-1); tindex=Ny-1:Ny-1:(Nx-1)*(Ny-1); % 境界条件により分かっている値を右辺に移項 U(lindex)=U(lindex)+betax*blvector; U(rindex)=U(rindex)+betax*brvector; U(bindex)=U(bindex)+betay*bbvector; U(tindex)=U(tindex)+betay*btvector; % 連立1次方程式を解いて差分解を求める U=AU\(AL\U(ap,:)); % 差分解を meshgrid に収める u(2:Ny,2:Nx)=reshape(U,Ny-1,Nx-1); % 境界条件を用いて境界上の格子点での値を求める u(:,1)=b_l(Y); u(:,Nx+1)=b_r(Y); u(1,:)=b_b(X); u(Ny+1,:)=b_t(X); % 解のグラフ figure; meshc(x,y,u); axis([a b c d -2 2]); title('graph of the approximate solution'); drawnow; shg; disp('何かキーを押して下さい。'); % 誤差を調べる 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 -2 2]); title('graph of the exact solution'); drawnow; shg end |
入手するには、ターミナルで
curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2d_nh.m |
cp -p poisson2d_nh.m ~/Documents/MATLAB |
>> poisson2d_nh(分割数はデフォールトの , が採用される。) |
>> poisson2d_nh(100,100) |
誤差は次のようになった。 となっていると考えている。
>> poisson2d_nh Nx=50, Ny=50, hx=0.020000, hy=0.020000, betax=2500.000000, betay=2500.000000 Nx=50, Ny=50, ||error||=5.924640e-05, ||u||=1.000000e+00 >> poisson2d_nh(100,100) Nx=100, Ny=100, hx=0.010000, hy=0.010000, betax=10000.000000, betay=10000.000000 Nx=100, Ny=100, ||error||=1.482114e-05, ||u||=1.000000e+00 >> poisson2d_nh(200,200) Nx=200, Ny=200, hx=0.005000, hy=0.005000, betax=40000.000000, betay=40000.000000 Nx=200, Ny=200, ||error||=3.705884e-06, ||u||=1.000000e+00 >> poisson2d_nh(400,400) Nx=400, Ny=400, hx=0.002500, hy=0.002500, betax=160000.000000, betay=160000.0000 00 Nx=400, Ny=400, ||error||=9.265347e-07, ||u||=1.000000e+00 |