数値積分とは面積を求める事です。
基本は台形則にあります。
求める面積を、多数の台形に分割します。
台形は計算できますので、合計すれば面積は求まります。
この時の多数とはどれくらいかがポイントです。
無限と言えば正解ですが、現実的には計算できません。
はじめに許容誤差を決めておきます。
分割する台形の数を徐々に増やします。
繰り返し、面積を求めて、１回前と比較します。
もしその差が、設定した許容誤差よりも小さくなれば、目標の近似になったと判断して終了します。
無限ではなくどこかで、終了する。これが基本です。
数値積分とは、ｘ軸に沿って複雑な形状をした「関数曲線」との間の面積を計算することです。
その方法は、「台形に近似する」で、これのｘ軸の幅を極限まで小さくすると正しい値になります。
極限までは、無理ですので、幅を狭くしてもほとんど変わらなくなれば十分に近似できたと考えます。
def fnf(x):=sin(x) 注：例
s:=0
repeat
h:=(b-a)/n
x:=a
f1:=abs(fnf(x))
i:=1
repeat
x:=a+h*i
f2:=abs(fnf(x))
s:=s+0.5*h*(f1+f2)
f1:=f2
i:=i+1
until i>=n
until ((a>b)+(n<=1))<=0
２関数積分も基本的な考えは同じです。
入れ子構造になります。
１関数と同じ台形則が、一番平凡で普通です。
f$:="(p^2*cos(q))" ２変数関数
input a,b
input c,d
input e,l
x1:=(b-a)*2
y1:=(d-c)*2
i:=1
repeat
x1:=x1*0.5
y1:=y1*0.5
y:=c
v:=0
i1:=2^(i-1)
j:=1
repeat
p:=a
q:=y
v2:=val(f$)
p:=b
v2:=v2+val(f$)
if i1:=1 then label1
x:=a
k:=1
repeat
c:=x+x1
p:=x
v2:=v2+2*val(f$)
k:=k+1
until k>=i1-1
label1
if j=1 or j=i1+1 then v2:=v2*0.5
v:=v+2*v2
y:=y+y1
j:=j+1
until j>=i1+1
v:=v*(b-a)*(d-c)/4^i
if i=1 then label2
if abs(v-v1)
v1:=v
i:=i+1
until i>=l
i:=i-1
label3
end
def fnf(x):=sin(x) */exanple/*
input a,b,n,e
s:=0
repeat
d:=b-a
x:=d/2+a
n1:=1
t:=abs(fnf(x)*d*2/3
s:=t+(abs(fnf(a))+abs(fnf(b)))*d/6
i:=2
repeat
s1:s
s:=(s-t/2)/2
n1:=n1*2)
c:=d/n1/2+a
t:=abs(fnf(c))
c1:=c
c2:=d/n1
j:=1
repeat
c1:=c1+c2
t:=t+abs(fnf(c1))
until j
s:=s+t
if abs(s-b1)
label1
until ((a>b)+(n<=1))>0
end
def fnf(x):=sin(x) ;example
input a,b,n,e
s:=0
if ((a>b)+(n<=1))>0 then l1
d:=b-a
x:=d/3+a
y:=d*2/3+a
n1:=1
t:=(abs(fnf(x))+abs(fnf(b)))+d/8
for i:=2 to n
s1:=s
s:=(s-t/3)/3
n1:=n1*3
c5:=d/n1/3*2+a
c4:=d/n1/3+a
t:=abs(fnf(c4))+abs(fnf(c5))
c6:=c4
c7:=c5
c2:=d/n1
for j:=1 to n1-1
c6:=c6+c2
c7:=c7+c2
t:=t+abs(fnf(c6))+abs(fnf(c7))
next t
t:=t*d*3/(n1*8)
s:=s+t
if abs(s-s1)
i:=i-1
label l2
end
label l1
end
def fnf(x):=sin(x) <---exsample-->
dim m(27);
data 1,0,2,-0.577350269,0.577350269,3,-0.707106781;
data 0,0.707106781,4,-0.794654472,-0.187592474,0.187592474;
data 0.794654472,5,-0.832497487,-0.374541109,0;
data 0.374541409,0.832497487,6,-0.866246818,-0.422518654;
data -0.266635402,0.266635402,0.422518654,0.866246818;
for i:=1 to 27;
read m(i);
next i;
label l1;
input a,b,e;
if a>=b then label l1;
a1:=a+b;
a2:=b-a;
i1:=1;
x:=a1*0.5;
s1:=2*abs(fnf(x));
for i:=2 to 6;
i1:=i1+i;
s:=0;
print i-1,s1;
for j:=1 to i;
j1:=i1+j;
x:=(a1+a2*m(j1))*0.5;
s:=s+abs(fnf(x));
next j;
s:=s*a2/m(i1);
if abs(s-s1)
next i;
i:=i+1;
label l2;
print "resalt=";s;spc$(3);"step=";i;
end
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))
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
f$:="p^2*cos(q)"
input a,b,c.d;
input e,l;
x1:=b-a;
y1:=d-c;
p:=a;
q:=c;
v1:=val(f$):
p:=b;
v1:=v1+val(f$);
q:=d;
v1:=v1+val(f$);
p:=a;
v1:=v1+val(f$);
v:=v1;
for i:=1 to l-1;
y2:=y1;
x2:=x1;
y1:=y1*0.5;
x1:=x1*0.5;
y:=c;
for j:=1 to 2^(i-1)+1;
q:=y;
x:=a+x1;
v2:=0;
for k:=1 to 2^(i-1);
p:=x;
v2:=v2+val(f$);
x:=x2;
next k;
if j=1 or j=2^(i-1)+1 then v2:=v2*0.5;
v:=v+v2*4;
y:=y+y2;
next j;
y:=c+y1;
for j:=1 to 2^(i-1);
q:=y;
x:=a;
v2:=0;
for k:=1 to 2^i+1;
p:=x;
if k=1 or k=2^I+1 then v2:=v2+2*val(f$);
goto l1;
v2:=v2+4*val(f$);
label l1;
x:=x+x1;
next k;
v:=v+v2;
y:=y+y2;
next j;
v3:=v*(b-a)*(d-c)/4^(i+1);
if abs(v3-v1)
next i;
i:=i-1;
l2;
print v3;
end;
def fnf(x):=sin(x);
dim m(27);
data 1,0,2,-0.577350269,0.577350269,3,-707106781;
data 0,0.707106781,4,-0.794654472,-0.187592474,0.187592474;
data 0.794654472,5,-0.832497487,-0.374541409,0;
data 0.374541409,0.832497487,6,-0.866246878,-0.422518654;
data -0.266635402,0.266635402,0.422518654,0.866246818;
for i:=1 to 27;
read m(i);
next i;
label l1;
input a,b,e,n,l;
if a>=b l1;
n1:=0;
for i:=1 to n;
n1:=n1+i;
next i;
dim n(n);
for i:=1 to n;
n(i):=m(n1+i);
next i;
a3:=(b-a)*2;
for i:=1 to l;
a5:=a;
a3:=a3*0.5;
s:=0;
for j:=1 to 2~(i-1);
a4:=a5;
a5:=a5+a3;
for k:=1 to n;
x:=(a4+a5+a3*n(k))*0.5;
s:=s+abs(fnf(x));
next k;
next j;
s:=s*a3/n;
if i:=1 then l2;
if abs(s-s1)
print i,d;
s1:=s;
next i;
i:=i-1;
labl l3;
print s,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;