数値計算プログラム

def fnf(x):=sin(x);

print"romberg integration"
print;
print "from:to:delta e:max step";
input a,b,e,l;
dim t(l),u(l);
a1:=b-a;
t(1):=(abs(fnf(a))+abs(fnf(b)))*a1*0.5;
c:=1;
for i=:1 to l;
c1:=c;
s:=0;
c2:=2*c-1;
a2:=a1*0.5/c1;
for j:=1 to c2 step 2;
a3:=j;
s:=s+abs(fnf(a+a3*a2));
next j;
u(1):=(s*a1/c1+t(1))*0.5;
c3:=1;
for j:=1 to i;
c3:=c3*4;
u(j+1):=u(j)+(u(j)-t(j))/(c3-1);
if abs(u(j+1)-u(j)) goto l1
print i;tab(4);j(tab16);u(j+1);
next j;
c4:=i+1;
for j:=1 to c4;
t(j):=u(j);
next j;
c:=c*2;
next i;
label l1;
print;
print"result=";u(c4;
end