「非同次 Neumann 境界条件下の熱方程式に対する差分法 -- MATLAB を使って数値計算」
heat2n_nh.m |
% heat2n_nh.m % 長方形領域における熱方程式 u_t=△u (非同次Neumann境界条件) を差分法で解く % 2024/8/24, written by mk. % 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2n_nh.m % 解説 https://m-katsurada.sakura.ne.jp/labo/text/heat2n-nonhomo.pdf function heat2n_nh function [A,B]=heat2n_mat0(Nx,Ny,lambdax,lambday,theta) 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; Kx=2*Ix-Jx; Ky=2*Iy-Jy; % column first A=kron(Ix,Iy)+theta*lambday*kron(Ix,Ky)+theta*lambdax*kron(Kx,Iy); B=kron(Ix,Iy)-(1-theta)*lambday*kron(Ix,Ky)-(1-theta)*lambdax*kron(Kx,Iy); % row first % A=kron(Iy,Ix)+theta*lambdax*kron(Iy,Kx)+theta*lambday*kron(Ky,Ix); % B=kron(Iy,Ix)-(1-theta)*lambdax*kron(Iy,Kx)-(1-theta)*lambday*kron(Ky,Ix); end % 長方形領域 (a,b)×(c,d) a=0; b=1; c=0; d=1; % Neumann境界値 b() 調和関数を選ぶとこれが平衡解 (定常解) function val = harmonic(x,y) val = (x-0.5).*(y-0.5); end % 初期値 (平衡解を境界値が0のもので摂動する) function val = iv(x,y) val = harmonic(x,y) + cos(pi*x).*cos(pi*y); end % 厳密解が分かる function val = exact(x,y,t) val = harmonic(x,y) + exp(- 2 * pi * pi * t) * cos(pi*x).*cos(pi*y); end % 長方形の下の辺での b_b() … -∂u/∂y(x,0,t) function val = b_b(x) val = 0.5 - x; end % 長方形の上の辺での b_b() … ∂u/∂y(x,1,t) function val = b_t(x) val = x - 0.5; end % 長方形の左の辺での b_l() … -∂u/∂x(0,y,t) function val = b_l(y) val = 0.5 - y; end % 長方形の右の辺での b_r() … ∂u/∂x(1,y,t) function val = b_r(y) val = y - 0.5; end % 誤差を測るためのノルム (最大値ノルム) function val = mynorm(a) [m,n]=size(a); val = norm(reshape(a,m*n,1),Inf); end % 格子への分割 Nx=50; Ny=50; hx=(b-a)/Nx; hy=(d-c)/Ny; theta=0.5; tau=0.5/(1/hx^2+1/hy^2); % 陽解法の場合のギリギリの刻み(もっと大きくてもOK) lambdax=tau/hx^2; lambday=tau/hy^2; display(sprintf("Nx=%d, Ny=%d, hx=%f, hy=%f, tau=%e, λx=%f, λy=%f", ... Nx, Ny, hx, hy, tau, lambdax, lambday)); % 差分方程式 A U^{n+1}=B U^n の行列 [A,B]=heat2n_mat0(Nx,Ny,lambdax,lambday,theta); % 対称にするためのスケーリング用の行列 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(full(er),inf)>=1 disp('Aが対称でない。') pause end % 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); % 初期値 u= iv(x,y); % 極限 ulimit = harmonic(x,y); % 初期値のグラフを描く disp('初期値のグラフ') meshc(x,y,u); axis([a b c d -4 4]); drawnow; disp('何かキーを押して下さい') pause % 境界値 左の辺、右の辺 blvector=b_l(Y(1:Ny+1)); brvector=b_r(Y(1:Ny+1)); % 境界値 下の辺、上の辺 bindex=1:Ny+1:(Nx+1)*(Ny+1); tindex=Ny+1:Ny+1:(Nx+1)*(Ny+1); bbvector=b_b(X(1:Nx+1)); btvector=b_t(X(1:Nx+1)); % どこまで計算するか Tmax=0.5; Nmax = ceil(Tmax/tau); % グラフを表示する時間の間隔 dt=0.005; skip=dt/tau; U=reshape(u,(Nx+1)*(Ny+1),1); disp('ループに入ります。何かキーを押してください。') t=0; maxerror=0; for n=0:Nmax % 注: これから計算するのが時刻tでの値 (最初が t=tau) U=B*U; % 移項 U(1:Ny+1)=U(1:Ny+1)+2*hx*lambdax*blvector; U(end-(Ny+1)+1:end)=U(end-(Ny+1)+1:end)+2*hx*lambdax*brvector; U(bindex)=U(bindex)+2*hy*lambday*bbvector; U(tindex)=U(tindex)+2*hy*lambday*btvector; U=S*U; % 連立一次方程式を解いて差分解を求める U=AU\(AL\U(ap,:)); t=(n+1)*tau; % 計算結果の表示 if mod(n,skip)==0 u=reshape(U,Ny+1,Nx+1); meshc(x,y,u); axis([a b c d -1 1]); drawnow; % 誤差の表示 error = mynorm(u-exact(x,y,t)); if error > maxerror maxerror = error; end display(sprintf("t=%f, ||error||=%e, ||u||=%e, ||u∞||=%e",... t, error, mynorm(u), mynorm(ulimit))); end end display(sprintf("max ||u(・,t)-u(・,∞)||=%e", mynorm(u-ulimit))); display(sprintf("Nx=%d,Ny=%d,max error=%e", Nx, Ny, maxerror)); end |
curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2n_nh.m cp -p heat2n_nh.m ~/Documents/MATLAB |
>> heat2n_nh >> heat2n_nh(100,100) |