微分は数値計算に向かないと言いながら、何故か最初はニュートン法から教科書は始まります。
習慣的に積分の前に微分をするのでしょうか。
ニュートン法というのは誤差が、終了設定値より小さくなるまで繰り返すだけです。
３乗根の解法
式は展開したものを使用します。(x*x*x=a)
d:=1e-11
x:=1
repeat
y:=x
x:=(x+3*a/(2*x*x+a/x))/2
until abs(x-y)>d
end
精度が悪い「数値微分」ですが、変化の緩い関数では可能です。
def fnf(x):=1/x example
x1 from
x2 to
x3 step
kizami
y1 from
y2 to
y3 step
x:=x1
repeat
y:=y1
f1:=(fnf(x+y)-fnf(x-y))/(2*y)
f2:=(fnf(x+y)-2*fnf(x-y)+fnf(x-y)/y^2
f:=fnf(x)
y:=y*y3
if y
label1
if x>x2 then end
until 0
input f$;
input h,n;
dim x(n),y(n);
input x(1),y(1);
for i=2 to n;
x:=x(i-1);
y:=y(i-1);
y(i):=y(i-1)+h*val(f$);
x(i):=x(i-1)+h;
next i;
label l1;
input x1,x2,x3;
if x1
h1:=x3/h;
h2:=(x1-x(1))/h+1:
h3:=(x2-x(1))/h+1;
for i=h2 to h3 step h1;
print x(1),y(i);
next i;
end;
input f$;
input h,n,m;
input x1,y1;
m1:=n/m;
x2:=x1;
y2:=y1;
for i:=1 to m1;
for j:=1 to m;
x:=x2;
y:=y2;
k1:=h*val(f$);
x:=x2+0.5*h;
y:=y2+0.5*k1;
k2:=h*val(f$);
y:=y2+0.5*k2;
k3:=h*val(f$);
x:=x2+h;
y:=y2+k3;
k4:=h*val(f$);
f:=(k1+2*(k2+k3)+k4)/6;
x2:=x2+h;
y2:=y2+f;
next j;
next i;
end
q1:=sqr(2);
input m;
input x1;
input y1;
m1:=n/m;
x2:=x1;
y2:=y1;
for i:=1 to m1;
for j:=1 to m;
x:=x2;
y:=y2;
p:=h*val(f$);
y:=y+0.5*p;
q:=p;
x:=x2+0.5*h;
p:=h*val(f$);
y:=y+(1-q1*0.5)*(p-q);
q:=(3/q1-2)*q+(2-q1)*p;
p:=h*val(f$);
y:=y+(1+1/q1)*(p-q);
q:=(2+q1)*p-(2+3/q1)*q;
x:=x2+h;
p:=h*val(f$);
y2:=y+p/6-q/3;
x2:=x2+h;
next j;
next i;
end;
input F$;
input h,e;
input n,l;
dim x(n),y(n),f(n);
input x(1),y(1);
for i:=1 to 5;
x:=x(i);
y:=y(i);
c1:=h*val(F$);
x:=x(i)+h/3;
y:=y(i)+c1/3;
c2:=h*val(f$);
x:=x(i)+h*0.4;
y:=y(i)+(6*c2+4*c1)/25;
c3:=h*val(f$);
x:=x(i)+h;
y:=y(i)+(15*c3-12*c2+c1)/4;
c4:=h*val(f$);
x:=x(i)+2*h/3;
y:=y(i)+(8*c4-50*c3+90*c2+6*c1)/81;
c5:=h*val(f$);
x:=x(i)+0.8*h;
y:=y(i)+(8*c4+10*c3+36*c2+6*c1)/75;
c6:=h*val(f$);
x(i+1):=x(i)+h;
y(i+1):=y(i)+(23*c1+125*c3-81*c5+125*c6)/192;
x:=x(i+1);
y:=y(i+1);
f(i+1):=val(f$);
next i;
for i:=6 to n-1;
p1:=1901*f(i)-2774*f(i-1)+2616*f(i-2)-1274*f(i-3)+251*f(i-4);
p:=y(i)+(h/720)*p1;
x(i+1):=x(i)+h;
for j:=1 to l;
x:=x(i+1);
y:=p;
f(i+1):=val(f$);
q1:=251*f(i+1)+646*f(i)-264*f(i-1)+106*f(i-2)-19*f(i-3);
q:=y(i)+h/720*q1;
if abs(p-q)
next j;
n:=i;
goto l2;
label l1;
y(i+1):=q;
x:=x(i+1);
y:=y(i+1);
f(i+1):=val(f$);
next i;
label l2;
input i1,i2;
if (i-1)*h
print x(i),y(i);
next i;
end;
k:=0;
j:=0;
input f$,h,e,n,l;
dim x(n+1),y(n+1),s(l+1),t(l+1);
input x(1),y(1);
for i:=1 to n;
x:=x(i);
y:=y(i);
t(i):=y+h*val(f$);
a:=2;
for j:=1 to l;
a1:=a;
x:=x(i);
y:=y(i);
for k:=1 to a;
f:=val(f$);
x:=x+h/a1;
y:=y+h/a1*f;
next k;
s(1):=y;
f1:=1;
for k:=1 to j;
f1:=f1*2;
s(k+1):=s(k)+(s(k)-t(k))/(f1-1);
if abs(s(k+1)-s(k))
j1:=j+1;
for k:=1 to j1;
t(k):=s(k);
next k;
a:=a+2;
next j;
goto l2;
label l1;
x(i+1):=x(i)+h;
y(i+1):=s(k+1);
next i
label l2;
input h1;
for i:=1 to n step h1/h;
print x(i),y(i);
next i;
end;
def fnf(x):=sin(x);
p:=3.14159565
i1:=0;
i:=0;
w:=0
a:=0;
b:=0;
label li;
input w,n;
input x1,x2,x3;
if n/2=int(n/2) then a:=n-0.5:b:=n+0.5; else a:=n+0.5:b:=n-0.5;
for x:=x1 to x2 step x3;
f:=0;
for i=-a to a step 2;
i1:=i*p;
f:=fnf(x+i1/w)/i1^2+f;
next i;
for i:-b to b step 2;
i1:=i*p;
f:=f-fnf(x+i1/w)/i1^2;
next i;
next x;
input i2;
if i2=0 then l1;
end;
２元
－－－－－
k:=0:
j:=0;
input f$,g$;
input h,e;
input n;
input l;
dim x(n+1),y(n+1),z(n+1),u(l+1),v(l+1),s(l+1),t(l+1);
input x(1),y(1),z(1);
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;
a1:=a;
x:=x(i);
y:=y(i);
z:=z(i);
for k:=1 to a;
f:=val(f$);
g:=val(g$);
x:=x+h/a1;
y:=y+h/a1*f;
z:=z+h/a1*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 l1;
next k;
j1:=j+1;
for k:=1 to j1;
t(k):=s(k);
v(k):=u(k);
next k;
a:=a*2;
next j;
go l2;
label l1;
x(i+1):=x(i)+h;
y(i+1):=s(k+1);
z(i+1):=u(k+1);
next i;
label l2;
input h1,h2;
if h2=0 then l3;
for i:=1 to n step h1/h;
print x(i),y(i);
next i;
end;
label l3;
for i:=1 to n step h1/h;
print x(i_),tab(3),y(i),z(i);
next i;
end;
r:=sqr(2);
input f$;
G$:="Z";
input e,i2;
input x1,y1,x2.y2;
label l1;
input h;
n:=(x2-x1)/h;
if int(n)<>n then print"re-enter":go l1;
dim x(n+1),y(n+1),a(i2+1),b(i2+1);
i1:=1;
label l7;
input z1;
label l6;
x3:=x1;
y3:=y1;
z3:=z1;
x(1):=x1;
y(1):=y1;
for i:=1 to n;
x:=x3;
y:=y3;
z:=z3;
p1:=h*val(f$);
p2:=h*val(g$);
y:=y+0.5*p1;
z:=z+0.5*p2;
q1:=p1;
q2:=p2;
x:=x3+0.5*h;
p1:=h*val(f$);
p2:=h*val(g$);
y:=y+(1-r*0.5)*(p1-q1);
z:=z+(1-r*0.5)*(p2-q2);
q1:=(3/r-2)*q1+(2-r)*p1;
q2:=(3/r-2)*q2+(2-r)*p2;
p1:=h*val(f$);
p2:=h*val(g$);
y:=y+(1+1/r)*(p1-q1);
z:=z+(1+1/r)*(p2-q2);
q1:=(2+r)*p1-(2+3/r)*q1;
q2:=(2+r)*p2-(2+3/r)*q2;
x:=x3+h;
p1:=h*val(f$);
p2:=h*val(g$);
y3:=y+p1/6-q1/3;
y(i+1):=y3;
x3:=x3+h;
z3:=z+p2/6-q2/3;
x(i+1):=x3;
next i;
if abs(y3-y2)
if i1=1 then a(1):=z1:b(i):=y3:i1:=2:go l3;
for o:=1 to i1-1;
if b(i)=y3 then l3;
next i;
a(i1):=z1;
b(i1):=y3;
for i:=1 to i1-1;
for j:=i+1 to i1;
if b(i)>=b(j) then l4;
w1:=b(i);
b(i):=b(j);
b(j):=w1;
w1:=a(i);
a(i):=a(j);
a(j):=w1;
label l4;
next j;
next i;
s1:=0;
s2:=1;
for i:=1 to i1;
s3:=1;
for j:=2 to i1;
if j=i then s3:=s3*(y2-b(j));go l5;
s3:=s3*(b(i)-b(j));
label l5;
next j;
s1:=s1+a(i)/s3;
s2:=s2*(y2-b(i));
next i;
z1:=s1*s2;
i1:=i1+1;
for i:=1 to i1-1;
if a(i)=z1 then l3;
next i;
gosub l6:go l6;
label l3;
gosub l6:go l7;
label l6;
for i:=1 to i1-1;
print a(i),b(i);
next i;
return;
label l2;
input n1;
for i:=1 to n+1 step n1/h;
print x(i),y(i);
next i;
end;
input p$,q$,e$;
input x1,y1,x2,y2;
input h;
input e;
input l;
n:=(x2-x1)/h;
if int(n)<>n then n:=int(n)+1:h:=(x2-x1)/n:print "step=",h;
dim x(n+1),y(n+1),z(n+1),a(n+1),u(l+1),v(l+1),s(l+1),t(l+1);
f$:="-z*val(p$)-y*val(q$)":g$:="z";
x(1):=x1;
y(1):=y1;
z(1):=0;
gosub l1;
for i:=1 to n+1;
a(i):=y(i);
next i;
f$:=f$+"+val(e$)";
gosub l1;
for i:=1 to n+1;
y(i):=(y2-y(n+1))*a(i)/a(n+1)+y(i);
next i;
input h1;
for i:=1 to n+1 step h1/h;
print x(i),y(i);
next i;
end;
label l1;
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;
z:=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)-c(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;
j1:=j+1;
for k:=1 to j1;
t(k):=s(k);
v(k):=u(k);
next k;
a:=a*2;
next j;
print "no converge";
go l3;
kabel l2;
x(i+1):=x(i)+h;
y(i+1):=s(k+1);
z(i+1):=u(k+1);
next i;
label l3;
return;
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;
初期値問題
オイラー・ロンバーグ法
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;
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 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 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 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 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;