Created
January 14, 2020 06:55
-
-
Save monhime/2a18b9826dc82848e0f1429dea3bd0b7 to your computer and use it in GitHub Desktop.
差分法による電位分布計算
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
N = 300; %r方向の格子数 | |
M = 300; %z方向の格子数 | |
RMAX = 30; %r座標の最大値 | |
ZMAX = 30; %z座標の最大値 | |
Conv = 1e-6; %収束判定誤差 | |
dr = RMAX / N; %r方向の刻み幅 | |
dz = ZMAX / M; %z方向の刻み幅 | |
phi = zeros(N,M); %0埋め | |
flag = -10*ones(N,M); %非境界部分の-10埋め | |
for i = 1:N | |
for j = 1:M | |
%格子点(i,j)の座標(r,z) | |
r = i*dr; | |
z = j*dz; | |
%ディリクレ条件(1kV): 高圧電極 | |
if ((r<1 && z > 16) || r^2+(z-14)^2<9) | |
flag(i,j)=-3; | |
phi(i,j)=1.0; | |
%ディリクレ条件(0kV): 低圧電極、接地面、r遠方境界 | |
elseif ((r<1 && z < 5) || r^2+(z-7)^2<9 || j==1 || i==N) | |
flag(i,j)=-3; | |
%z方向ノイマン条件 | |
elseif (j == M) | |
flag(i,j)=-1; | |
%r方向ノイマン条件 | |
elseif (i == 1) | |
flag(i,j)=-2; | |
end | |
end | |
end | |
MaxPhi = 1.0; | |
MaxErr = 1e20; | |
loop = 0; | |
while MaxErr > Conv | |
Prev_phi=phi; | |
for i = 1:N | |
for j = 1:M | |
r = i*dr; %格子点(i,j)の座標(r,z) | |
%非境界のphi | |
if (flag(i,j)==-10) | |
phi(i,j)=0.25*(Prev_phi(i+1,j)+Prev_phi(i-1,j)+Prev_phi(i,j+1)+Prev_phi(i,j-1)+0.5*dr*(Prev_phi(i+1,j)-Prev_phi(i-1,j))/r); | |
%z方向ノイマン境界 | |
elseif (flag(i,j)==-1) | |
phi(i,j)=0.25*(Prev_phi(i+1,j)+Prev_phi(i-1,j)+2*Prev_phi(i,j-1)+0.5*dr*(Prev_phi(i+1,j)-Prev_phi(i-1,j))/r); | |
%r方向ノイマン条件 | |
elseif (flag(i,j)==-2) | |
phi(i,j)=0.25*(2*Prev_phi(i+1,j)+Prev_phi(i,j+1)+Prev_phi(i,j-1)); | |
end | |
end | |
end | |
%前回の行列との誤差。MaxPhiで規格化 | |
Err = (abs(phi-Prev_phi))/MaxPhi; | |
MaxErr=max(max(Err));%行列の最大の要素 | |
loop = loop + 1; | |
%1000ループごとにMaxMaxErr表示 | |
if rem(loop,1000)==0 | |
disp(loop); | |
disp(MaxErr); | |
end | |
end | |
rdata=linspace(dr,dr*N,N); | |
zdata=linspace(dz,dz*M,M); | |
createfigure(rdata,zdata,phi'); | |
function createfigure(xdata1,ydata1,zdata1) | |
figure1 = figure; % figure を作成 | |
axes1 = axes('Parent',figure1);% axes を作成 | |
hold(axes1,'on'); | |
% contour を作成 | |
[c1,h1] = contour(xdata1,ydata1,zdata1,'TextStep',0.1,'LevelStep',0.01); | |
clabel(c1,h1); | |
ylabel({'z方向'}); % ylabel を作成 | |
xlabel({'r方向'}); % xlabel を作成 | |
box(axes1,'on'); | |
axis(axes1,'tight'); | |
% 残りの座標軸プロパティの設定 | |
set(axes1,'BoxStyle','full','FontName','Times New Roman','FontSize',24,'Layer','top','XGrid','on','XMinorGrid','on','XTick',[5 10 15 20 25 30],'YGrid','on','YMinorGrid','on','YTick',[5 10 15 20 25 30]); | |
colorbar(axes1); % colorbar を作成 | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment