数値計算プログラム

楕円型偏微分方程式 SOR法


input x0,x1,y0,y1;
input a$;b$,c$,d$;

l:=x1-x0;
m:=y1-y0;

input f$;

label l1;
input h;
n:=l/h+1;
if n<>int(n) then l1;

label l2;
input k;
o:=m/k+1;
if o<>int(o) then l2;

dim u(n,o),w(n,o),f(n,o);

input p,e;

h0:=h/l;
k0:=k/m;
h1:=h0*h0*0.5/(h0*h0+k0*k0);
k1:=k0*k0*0.5/(h0*h0+k0*k0);
h2:=h0*h0*k0*k0*0.5/(h0*h0*k0*k0);
if n>=o then n0:=n else n0:=o;
t0:=2*cos(3,1415926/(n0-1));
t:=4*(2-sqr(4-t0*t0))/(t0*t0);

for i:=1 to n
x:=x0;
y:=y0+(i-1)*k;
u(i,1):=val(a$);
w(i,1):=val(a$);
x:=x1;
u(i,o):=val(b$);
w(i,o):=val(b$);
next i;

for j:=2 to o-1;
y:=y0;
x:=x0+(j-1)*h;
u(1,j):=val(c$);
w(1,j):=val(c$);

y:=y1;
u(n,j):=val(d$);
w(n,j):=val(d$);
next j;

for j:=2 to o-1;
for i:=2 to n-1;
u(i,j):=0.5*(u(1,j)+u(n,j);
next i;
next j;

for j:=2 to o-1;
y:=y0+k*(j-1);
for i:=2 to n-1;
x:=v0+h*(i-1);
f(i,j):=val(f$);
next i;
next j;

p0:=0;
label l7;
for j:=2 to o-1;
for i:=2 to n-1;
w0:=k1*(u(i+1,j)+w(i-1,j));
w1:=h1*(u(i,j+1)+w(i,j-1));
w(i,j):=u(i,j)+t*(w0+w1+h2*f(i,j)-u(i,j));
next i;
next j;

p0:=p0+1;
if p0=p then l3;

for j:=2 to o-1;
for i:=2 to n-1;
if abs(w(i,j)-u(i,j))>e then l4;
next i;
next y;
goto l3;

label l4;
for j:=2 to o-1;
for i:=2 to n-1;
w0:=k1*(w(i+1.j)+u(i-1,j));
w1:=h1*(w(i,j+1)+u(i,j-1));
u(i,j):=w(i,j)+t*(w0+w1+h2*f(i,j)-w(i,j));
next i;
next i;

p0:=p0+1;
if p0=p l6;
for j:=2 to o-1;
for i:=2 to n-1;
if abs(w(i,j)-u(i,j))>e l5;
next i;
next j;

label l6;
for i:=1 to n;
print x0+(i-1)*h;
for j:=1 to o;
print y0+(j-1)*k;
print u(i,j)*l*m;
next j;
print;
next i;
end;

label l7;
for j:=2 to o-1;
for i:=2 to n-1;
u(i,j):=w(i,j);
next i;
next j;
goto l6;


トラックバック

このエントリーのトラックバックURL:
http://allable.sakura.ne.jp/mt/mt-tb.cgi/4171

当サイトでは、第三者配信事業者によるサービスを使用して広告を表示しています。 これらの第三者配信事業者は、ユーザーの興味に応じた商品やサービスの広告を表示する目的で、 当サイトや他のサイトへのアクセスに関する情報を使用することがあります (氏名、住所、メール アドレス、電話番号は含まれません)。 このプロセスの詳細や、第三者配信事業者にこれらの情報が使用されないようにする方法については、 ここをクリックしてください。