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;
input l;
input g$;
input f$;
label l1;
input h;
n:=l/h+1;
if n<>int(n) then l1;
dim u(n,n),w(n,n),f(n,n);
input m,e;
k:=h/l;
for i:=1 to n;
u(1,i):=0;
u(n,i):=0;
u(i,n):=0;
w(1,i):=0;
w(n,i):=0;
w(i,n):=0;
x:=(i-1)*k*l;
u(i,1):=val(g$);
w(1,1):=val(g$);
next i;
t:=2*cos(3.1415926/(n-1);
w0:=(2-sqr(4-t*t))/(t*t);
k0:=k*k;
for i:=2 to n-1;
for j:=2 to n-1;
u(i,j):=0;
next j;
nexy i;
m1:=0;
for j:=1 to n;
y:=(j-1)*k;
for i:=1 to n;
x:=(i-1)*k;
f(i,j):=val(f$);
next i;
next j;
label l4;
for j:=2 to n-1;
for i:=2 to n-1;
w:=w(i-1,j)+u(i+1,j)+w(i,j-1)+u(i,j+1)-4*u(i,j)+k0*f(i,j);
w(i,j):=u(i,j)+w0*w;
next i;
next j;
m1:=m1+1;
if m1=m then l2;
for j:=2 to n-1;
for i:=2 to n-1;
if abs(u(i,j)-w(i,j))>e then l3;
next i;
next j;
l2;
label l3;
for j:=2 to n-1;
for i:=2 to n-1;
u(i,j):=w(i,j);
next i;
next j;
l4;
label l2;
for i:=1 to n;
print (i-1)+h;
for j:=1 to n;
print (j-1)*h;
print (i-1)*h,u(i,j)*l*l;
next j;
print:
next i;
end;
input h;
input r;
k:=r*h+h;
input m;
input n;
dim a(n),u(n,m);
input x0;
input a$;
label l1;
input x1;
if (x1-x0)<>h*(n-1) then print"error";goto l1;
input b$;
input t0;
input c$;
c1:=2+2/r;
c2:=2-2/r;
a(2):=c1;
n1:=n-1;
for i:=3 to n1;
a(i):=c1-1/a(i-1);
next i;
for i:=1 to n;
t:=t0;
x:=(i-1)*h+x0;
u(i-1):=val$(c$);
next i;
for j:=1 to m;
x:=x0;
t:=t0+(j-1)*k;
u(1,j):=val$(a$);
x:=x1;
u(n,j):=val(b$);
next j;
for j:=2 to m;
for i:=2 to n1;
u(i,j):=u(i-1,j-1)-c2*u(i,j-1)+u(i+1,j-1);
next i;
u(2,j):=u(2,j)+u(1,j);
for i:=3 to n1;
u(i,j):=u(i,j)+u(i-1,j)/a(i-1);
next i;
for i:=n1 to 2 step -1;
u(i,j):=(u(i,j)+u(i+1,j))/a(i);
next i;
next j;
for j:=1 to m;
print"T=";t0+(j-1)*k;
for i:=1 to n;
print" x=";x0+(i-1)*h;"u=";u(i,j);
next i;
print;
next j;
end;
input h,r,;
r1:=r/(1+r)*0.5;
k:=r+h*h;
input m,n;
dim u(n,m(,b(n-2),w(n-1);
input x0.a$,x1,b$,t0,c$,l,e;
for j:=1 to m;
x;=x0;
t:=t0+k*(j-1);
u(1,j):=val(a$);
x:=x1;
u(n,j):=val(b$);
next j;
t:=t0;
for i:=2 to n-1;
x:=x0+h*(i-1);
u(i,1):=va;l$(c$);
for j:=2 to m;
u(i,j):=0;
next j;
w(i-1):=0:
next i;
w(n-1):=u(n,m);
w0:=cos(3.15159265/(m-1)/(1+r)*r;
w1:=2/(1+sqr(1-w0*w0));
w2:=w1-1;
w3:=r1*w1;
for j:=2 to m;
for i:=2 to n-1;
b(i-1):=u(i-1,j-1)+2*(1/r-1)*u(i,j-1)+u(i+1,j-1);
b(i-1):=b(i-1)*r1;
next i;
l1:=0;
label l3;
for i:=2 to n-1;
u(i,j):=w3*(u(i-1,j)+w(i)+b)i-1))-w2*w(i-1);
next i;
l1:=l1+1;
if l1=l then l1;
for i:=2 to n-1;
if abs(u(i,j)-w(i-1))>e then l2;
next i;
goto l1;
label l2;
for i:=2 to n;
w(i-1):=u(i,j);
next i;
goto l3;
label l1;
for i:=1 to n-1;
w(i):=0;
next i;
next j;
print "result";
end;
input h,r,;
r1:=r/(1+r)*0.5;
k:=r+h*h;
input m,n;
dim u(n,m(,b(n-2),w(n-1);
input x0.a$,x1,b$,t0,c$,l,e;
for j:=1 to m;
x;=x0;
t:=t0+k*(j-1);
u(1,j):=val(a$);
x:=x1;
u(n,j):=val(b$);
next j;
t:=t0;
for i:=2 to n-1;
x:=x0+h*(i-1);
u(i,1):=va;l$(c$);
for j:=2 to m;
u(i,j):=0;
next j;
w(i-1):=0:
next i;
w(n-1):=u(n,m);
for j:=2 to m;
for i:=2 to n-1;
b(i-1):=u(i-1,j-1)+2*(1/r-1)*u(i,j-1)+u(i+1,j-1);
b(i-1):=b(i-1)*r1;
next i;
l1:=0;
label l3;
for i:=2 to n-1;
u(i,j):=(u(i-1,j)+w(i))*r1+b(i-1);
next i;
l1:=l1+1;
if l1=l then l1;
for i:=2 to n-1;
if abs(u(i,j)-w(i-1))>e then l2;
next i;
goto l1;
label l2;
for i:=2 to n;
w(i-1):=u(i,j);
next i;
goto l3;
label l1;
for i:=1 to n-1;
w(i):=0;
next i;
next j;
print "result";
end;
円周率の計算例
input n;
a:=0;
b:=Rnd(7);
for i:=1 to n;
b:=(Rnd(0)~2+Rnd(0)~1);
if Sqr(b)<1 then a:=a+1;
nwxt i;
a:=a*4/n;
print a;
ens;
input a$,b$,c$,d$,x1,y1,x2,y2,h;
n:=(x2-x1)/h;
if int(n)=0 then l1;
n:=int(n)+1;
h:=(x2-x1)/n;
print h;
label l1;
input e,l;
dim x(n+1),y(n+1),z(n+1),u(l+1),v(l+1),s(l+1),t(l+1),a(10),b(56);
F$:="(val(d$)-val(b$)*z-val(c$))/val(a$)";
g$:="z";
x(1):=x1;
y(1):=y1;
z(1):=0;
i1:=1;
label l5;
for i:=1 to n;
x:=x(i);
y:=y(i);
z:=z(i);
t(1):=y+h*val(f$);
v(1):=z+h*val(g$);
a:=2;
for j:=1 to l;
x:=x(i);
y:=y(i);
z:=z(i);
for k:=1 to a;
f:=val(f$);
g:=val(g$);
x:=x+h/a;
y:=y+h/a*f;
z:=z+h/a*g;
next k;
s(1):=y;
u)1):=z;
f1:=1;
for k:=1 to j;
f1:=f1*2;
s(k+1):=s(k)+(s(k)-t(k))/(f1-1);
u(k+1):=u(k)+(u(k)-v(k))/(f1-1)
if abs(s(k+1)-s(k))<=s(k)*e and abs(u(k+1)-u(k))<=u(k)*e then l2;
next k;
for k:=1 to j+1;
t(k):=s(k);
v(k):=u(k);
next k;
a:=a*2;
next j;
goto l3;
label l2;
x(i+1):=x(i)+h;
y(i+1):=s(k+1);
z(i+1):=u(k+1);
next i;
label l3;
if abs(y(n+1)-y(2))
if i1=1 then i1:=2:z(i):=2*y2-y(n+1):a(1):=y(n+1):b(1):=0::goto l5;
for i:=1 to i1-1;
if a(i)=y(n+1) then label l6;
next i;
a(i1):=y(n+1);
b(i1*(i1-1)/2+1):=x(1);
for i:=2 to i1;
i0:=i1*(i1-1)/2;
i2:=i-1;
b(i0+i):=((a(i1)-y2)*b(i2*i2+1)/2)-(a(i2)-y2)*b(i0+i2))/(a(i1)-a(i2));
next i;
i1:=i1+1;
i3:=b(i1*(i1-1)/2);
for i:=1 to i1-2;
if b(i*(i+1)/2)=i3 then l6;
next i;
x(1):=i3;
goto l5;
label l6;
for i:=1 to i1-1;
print b(i*(i+1)/2),a(i);
next i;
input z(1);
goto l5;
input h1;
h2:=int(h1/h);
if h2<1 then h2:=1;
for i:=1 to n+1 step h2;
print x(i),y(i);
nexy i;
end;
初期値問題
オイラー・ロンバーグ法
input a$;
input b$;
input c$;
input d$;
input x0,y0,z0,v0,w0;
input h;
input n;
input e;
input l;
if n>254 then n:=254;
if l>255 then l:=254;
dim y(n+1),z(n+1),v(n+1),w(n+1,a(l+1),b(l+1),c(l+1),d(l+1),o(l+1),p(l+1),q(l+1),r(l+1;
x1:=x0;
y(1):=y0;
z(1):=z0;
v(1):=v0;
w(1):=w0;
for i:=1 to n;
x:=x1;
y:=y(i);
z:=z(i);
v:=(vi);
w:w(i);
a1:=2;
a(1):=y+h*val(z$);
b(1):=z+h*val(b$);
c(1):=v+h*val(c$);
d(1):=w+h*val(d$);
for j:=1 to l;
x:=x1;
y:=y(i);
z:=z(i);
v:=v(i);
w:=w(i);
for k:=1 to a1;
a:=val(a$);
b:=val(b$);
c:=val(c$);
d:=val(d$);
x:=x+h/a1;
y:=y+h/a1*a;
z:=z+h/a1*b;
v:=c+h/a1*c;
w:=w+h/a1*d;
nezt k;
o(1):=y;
p(1):=z;
q(1):=v;
r(1):=w;
a2:=1;
for k:=1 to j;
a2:=a2*2;
a3:=a2-1;
o(k+1):=o(k)+(o(k)-a(k))/a3;
o1:=abs(o(k+1)-o(k));
p(k+1):=p(k)+(p(k)-b(k))/a3;
p1:=abs(p(k+1)-p(k));
q(k+1):=(q(k)+(q)k)-c(k)/a3;
q1:=abs(q(k+1)-q(k));
r(k+1):=r(k)+(r(k)-d(k))/a3;
r1:=abs(r(k+1)-r(k));
if o1<=o(k)*e and p1<=p(k)*e and q1<=q(k)*e and r1<=r(k)*e then l1;
next k;
for k:=1 to j+1;
a(k):=o(k);
b(k):=p(k);
c(k):=q(k);
d(k):=r(k);
next k;
a1:=a1*2;
next j;
go l2;
label l1;
x1=x1+h;
y(i+1):=o(k+1);
z(i+1):=p(k+1);
v(i+1):=q(k+1);
w(i+1):=r(k+1);
next i;
label l2;
print x,y(i),z(i),z(i),w(i);
end;
F$:="x^2*cos(y)"
k:=0;
j:=0:
input a,b,c,d;
input e,l;
dom t(l),u(l);
x1:=b-a;
y1:=d-c;
x:=a;
y:=c;
t(1):=val(f$);
x:=b;
t(1):t(1)+val(f$);
y:d;
t(1):=t(1)+val(f$);
x:=a;
t(1):=t(1)+val(f$);
h:＝１；ｖ：＝ｔ（１）；
for i:=1 to l-1;
h1:=h;
h2:=2*h-1;
y2:=y1;
x2:=x1;
y1:=y1*0.5;
x1:=x1*0.5;
y:=c-y2;
i1:=2^(i-1);
for j:=1 to i1+1;
x:=a+x1;
v2:=0:
y:=y+y2;
for k:=1 to i1;
v2:=v2+val(f$);
x:=x+x2;
next k:
if j=1 or j=i1+1 then v:=v+v2*2;:go l1
v:=v+v2*4;
label l1;
nwxt j;
y:=c+y1;
for j:=1 to i1;
x:=a-x1;
v2:=0;
i2:=i1*2+1;
for k:=1 to i2;
x:=x+x1;
if k=1 or k=i2 then v2:=v2+2*val(f$);:go l2;
v2:=v2+4*val(f$);
label l2;
next k;
v:=v+v2;
y:=y+y2;
next j;
u(1):=v*(b-a)*(d-c)/4^(i+1);
h3:=1;
for j:=1 to i;
h3:=h3*4;
u(j+1):=u(j)+u(j)-t(j))/(h3-1);
if abs(u(j+1)-u(j))
next j;
j:=j+1;
h4:=i+1;
for k:=1 to h4;
t(k):=u(k);
next k;
h:=h*2;
next i;
label l3;
print u(h4),spc$(3),i+1,j;
end;
input p$;
input q$;
input F$;
input x1,y2,x2,y2;
input h;
n:=(x2-x1)/h
if int(n)=n then l1;
n:=int(n)+1;
h:=(x2-x1)/n;
print h;
label l1;
n:=n+1;
dim a(n,3),b(n),y(n),s(n),r(n);
x:=x1;
r(1):=0;
r(n):=0;
for i:=2 to n-1
v:=x+h;
a(i,1):=1-h*val(p$)/2;
a(i,2):=-(2-h*h*val(q$));
a(i,3):=1+h*val(p$)/2;
b(i):=h*h*valf$);
r(i):=0;
bext i;
x:=x1;
a(1,1):=0;
a(1,2):=-(2-h*h*val(q$));
a(1,3):=1+h*val(p$)/2;
b(1):=h*h*val(f$)-(1-h*val(p$)/2)*y1;
x:=x2;
a(n,1):=1-h*val(p$)/2;
a(n,2):=-(2-h*h*val(q$));
a(n,3):=0;
b(n):=h*h*val(f$)-(1+h*val(p$)/2)*y2;
s(1):=0;
for i:=2 to n;
s(i):=b(i-1);
for j:=i-2 to i-1;
if j<1 or j>n then l2;
s(i):=s(i)-s(j)*a(i-1,j-i+3);
label l2;
next j;
s(i):=s(i)/a(i-1,3);
next i;
y(1):=1;
for i:=2 to n;
y(i):=b(i-1);
for j:=i-2 to i-1;
if j<1 or j>n then l3;
y(i):=y(i)-y(j)*a(i-1,j-i+3);
label l3;
next j;
y(i):=y(i)/a(i-1,3);
r(i):=y(i)-s(i);
next i;
w:=b(n)-a(n,1)*s(n-1)-a(n,2)*s(n);
y(1):=w/(a(n,1)*r(n-1)+a(n,2)*r(n));
for i:=2 to n;
y(i):=y(1)*r(i)+s(i);
next i;
end;