「数値解析」は数学かどうかは意見が分かれるでしょう。
ここでは一応数学とします。
それは、デジタル値を相手にする所からすでに道が分かれますが、それでいて無限という概念がかならずついて回ります。
限りなく無限に近い、いいかえれば有限で計算を終える事は、近似を意味します。
はじめから、近似を前提にした数学。
非常に微妙な位置にあると言えます。
数値計算は近似です。
近似計算には必ず誤差があります。
しかし、誤差の大きさは計算対象で大きく異なります。
微分は、式の上では必ず計算できますが、線の傾きですのでわずかなずれで大きな誤差が発生します。
積分は、式では計算出来ないものも多くあります。しかし、元々は面積計算ですから、目分量で計算できる形に変形しても意外と誤差は少ないです。
数値積分は数値計算向きです。
補間法は、他の解法の部分処理として使用される事が多いと思っています。
逆にいえば、基本的な手法です。
dim x(n),y(n)
s:=0
f:=1
i:=1
repeat
d:=1
j:=1
repeat
if j=i then d:=d*(x-x(j)):go label1
d:=d*(x(i)-x(j)
label1
until j>=n
s:=s+y(i)/d
f:=f*(x-x(i))
until i>=n
end
ソーティングは数値解析ではないが、有用性は同等か上回ります。
シェル・ソートは、やや変則的な方法ですが、かなり良い結果がでます。
dim x(n)
n0:=n
label3
repeat
no:=int(no/2)
label2
n1:=0
n2:=n-n0
i:=1
repeat
n3:=i+n0
if x(i)<=x(n3) then label1
n4:=x(i)
x(i):=x(n3)
x(n3):=n4
n1:=1
label1
until i>=n2
if n1=1 label2
go label3
end
記憶容量は必要ですが、ソート速度が速いほうほうです。
分割してソートしてから、まとめてゆく方法です。
dim s(n,2)
n0:=1
s(1,1:=1
s(1,2):=n
label6
n1:=s(n0,1)
n2:=s(n0,2)
n0:=n0-1
label5
n3:=n1
n4:=n2
n5:=x((n1+n2)/2)
label1
if x(n3)
label2
if n5
if n4-n1<=n2-n3 then label3
if n3
label3
if n1
label4
if n1
return
dim r(n),j(n);
for i:=1 to n;
r(i):=0:
j(i):=0:
next i;
for i:=1 to 10
r(i):=20:
next i:
gosub l1:
for i:=1 to 40;
r(i):=sqr(r(i)*r(i)+j(i)*j(i));
next i;
j:=0;
k:=0;
for i:=1 to 40;
if r(i)>j then j:=r(i);
if r(i)
for i:=1 to 40;
print i,spc$(int(0.5+25*(r(i)-k)/(j-k)));
next i;
end
label l1;
n1:=1;
n2:=n;
p:=6.28318/n;
n2:=n2/2;
k1:=0;
label l2;
x:=0;
k2:=k1+n2-1;
for i:=k1 to k2;
h:=i+n2;
v:=cos(p*x);
w:=sin(p*x);
if f=1 then w:=-w;
r1:=r(i+1);
j1:=j(i+1);
r(i+1):=r1+r(h+1);
j(i+1):=j1+j(h+1);
q:=(r1-r(h+1))*v+(j1-j(h+1))*w;
j(h+1):=(j1-j(h+1))*v-(r1-r(h+1))*w;
r(h+1):=q;
x:=x+n1;
next i;
k1:=k1+n2+n2;
if k1
n1:=n1+n1;
goto l4;
label l3;
for i:=0 to n-1;
n2:=1;
h:=0;
n3:=n*0.5;
x:=1;
label l7;
n1:=x-n3;
if n1<0 then l5;
h:=h+n2;
x:=n1;
label l5;
if n3=1 then l6;
n3:=n3/2;
n2:=n2+n2;
goto l7;
label l6;
if i>=h then l8;
r1:=r(h+1);
j1:=j(h+1);
r(h+1):=r(i+1);
j(h+1):=j(i+1);
r(i+1):=r1;
j(i+1):=j1;
label l8;
next i;
return;
Pk,n(X)=((Xk-X)Pn-1,n-1-(Xn-1-X)Pk,n-1)/(Xn-Xn-1)
dim x(22),y(255);
labei li;
if h>22 then print"too many datas" go l1;
for i:=1 to k;
input x(i);
input y(i*(i-1)/2+1);
next i;
i1:=2;
for i:=i1 to k;
for j:=2 to i;
i0:=i*(i-1)/2;
j0:=j-1;
y(i0+j)=((x(i)-x)*y(j0*(j0+1)/2)-(x(j0)-x)*y(i0+j0)))/(x(i)-x(j0));
next j;
next i;
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;