連立方程式は行列の解になります。
一番オーソドックスなのが、「ガウスのはきだし法」です。
dim a(n,n+1) 不定の処理は省いています
i:=1
j:=1
repeat i
repeat j
input a(i,j)
until j>=n+1
until i>=n
k:=1
repeat k
w:=a(k,k)
j:=k
repeat j
a(k,j)*=a(k,j)/w
until j>=n+1
i:=1
repeat i
w:=a(i,k)
j:=k
repeat j
a(i,j):=a(i,j)-w*a(k,j)
until j>=n+1
until i>=n or i:=k
until k>=n or w
input n;
dim a(n,n+1);
for i:=1 to n;
for j:=1 to n+1;
input a(i,j);
next j;
next i;
for k:=1 to n;
label l5;
w:=a(k,k);
if w:=0 then l1;
for i:=1 to n;
if i=k then l2;
t:=a(i,k);
for j:=1 to n+1;
a(i,j):=w*a(i,j)-t*a(k,j);
next j;
label l2;
next i;
next k;
for i:=1 to n;
print a(i,n+1)/a(i/i);
next i;
end;
label l1;
if k=n then l3;
for h:=k+1 to n;
if a(h,h)<>0 then print h:goto l4
next h;
goto l3;
label l4;
for l:=1 to n+1;
v:=a(k,l);
a(k,l)=a(h,l);
a(h,l):=v;
next l;
goto l5;
label l3;
if a(k,n+1)=0 then print "不定";
end;
input n,e.j;
l1=0;
dim a(n,n+1),b(n),c(n);
for i:=1 to n;
c(i):=o;
for j:=1 to n+1;
input a(i,j);
next j;
next i;
for i:=1 to n;
label l8;
w:=a(i,i);
if w:=0 then l1;
for j:=1 to n+1;
a(i,j):=a(i,j)/w;
next j;
next i;
label l5;
for i:=1 to n;
b(i):=0;
for j:=1 to n;
if i=j then l2;
b(i):=b(i)-c(j)*a(i,j);
label l2;
next j;
b(i):=b(i)+a(i,j);
next i;
l1:=l1+1;
if l1=l then l3;
next i;
for i:=1 to n;
if abs(b(i)-c(i))>e then l4;
next i;
label l3;
for i:=1 to n;
print b(i);
next i;
end;
label l4;
for j:=1 to n;
c(j):=b(j);
next j;
goto l5;
label l1;
if i=n then l6;
for k:=i+1 to n;
if a(h,h)<>0 then i,h:goto l7;
next h;:goto l6;
label l7;
for l:=1 to n;
v:=a(i,l);
a(i,l):=a(h,l);
a(h,l):=v;
next l;
goto l8;
label l7;
print "funou";
end;
input n;
dim a(n,n+1),x(n);
input e;
for i:=1 to n;
for j:=1 to n;
input a(i,j);
next j;
next i;
for i:=1 to n;
input a(i,i+1);
next i;
n1:=n+1;
n2:=n-1;
for k:=1 to n2;
a1:=a(k,k);
for j:=k to n1;
a(j,j):=a(k,j)/a1;
next j;
for i:=k+1 to n;
a2:=a(i,k);
for j:=k to n1;
a(i,j):=a(i,j)-a2*a(k,j);
next j;
next i;
next k;
if abs(a(n,n1))<=e then l1;
a3:=a(n,n);
a(n,n)=a(n,n)/a3;
a(n,n1):=a(n,n1)/a3;
x(n):=a(n,n1);
for i:=1 to n2;
k:=n-i;
x(k):=a(k,n1);
for j:=k to n2;
x(k):=x(k)-a(k,j+1)*x(j+1);
next j;
next i;
for i:=1 to n;
print x(i);
next i;
end;
label l1;
print "error";
end;
inout n;
e:=1e-6
dim a(n,n),l(n,n),u(n,n),b(n),y(n);
for i:=1 to n;
for l:=1 to n;
imput a(i,j);
l(i,j):=0;
u(i,j):=0;
next j;
next i;
for i:=1 to n;
input b(i);
y(i):=0:
next i;
for j:=1 to n;
for i:=1 to n;
if i
for k:=1 to j;
w:=w+l(i,k)*u(k,j);
next k;
l(i,j):=a(i,j)-w;
goto l2;
label l1;
w=0;
for k:=1 to j-1;
w:=w+l(i,k)*u(k,j);
next k;
u(i,j):=(a(i,j)-w)/l(i,i);
label l2;
next i;
next j;
for i:=1 to n;
for j:=1 to i;
if i=j then l3;
b(i):=b(i)-y(j)*l(i,j);
label l3;
next j;
y(i):=b(i)/l(i,i)
next i;
for i:=n to 1 step -1;
for j:=i to n;
if i:=j then l4;
y(i):=y(i)-y(i)*u(i,i);
label l4;
next j;
next i;
for i:=1 to n;
print y(i);
next i;
end;
行列の対角線の３項のみ「０」以外の方程式。
input n;
dim a(n,3),b(n),x(n),s(n),r(n);
for i:=1 to n;
for j:=i-1 to i+1;
if j<1 or j>n then l1;
input a(i,j-i+2);
label l1;
next j;
next i;
for i:=1 to n;
input b(i);
r(i):=0;
next i;
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;
x(1):=1;
for i:=2 to n;
x(i):=b(i-1);
for j:=i-2 to i-1;
if j<1 or j>n then l3;
x(i):=x(i)-x(j)*a(i-1,j-i+3);
next j;
x(i):=x(i)/a(i-1,3);
r(i):=x(i)-s(i);
next i;
x1:=b(n)-a(n,1)*s(n-1)-a(n,2)*s(n);
x(1):=x1/(a(n,1)*r(n-1)+a(n,2)*r(n));
for i:=2 to n;
x(i):=x(1)*r(i)+s(i);
next i;
end;
def fnf(x):=x*x-exp(-x);
input l,e;
label l1;
input a,b;
if fnf(a)*fnf(b)>=0 then l1;
for i:=1 to l;
a1:=fnf(a);
b1:=fnf(b);
c:=(a*b1-b*a1)/(b1-a1);
f:=fnf(c);
if abs(f)
a:=c;
label l3;
next i;
i:=t-1;
label l2;
print c,f,i;
end;
注：関数は例。
input n;
dim a(n,n),x(n),b(n);
input e;
for i:=1 to n;
for j:=1 to n;
input a(i,j);
next j;
next i;
for i:=1 to n;
input b(i);
x(i):=i;
next i;
for k:=1 to n;
if k<>n then l1;
if abs(a(k,k)
label l1;
i1:=k;
a1:=abs(a(k,k));
k1:=k+1;
for i:=k1 to n:
a2:=abs(a(i,k));
if a2<=a1 then l4;
i1:=i;
a1:=a2;
label l4;
next i;
i2:=k;
a3:=abs(a(k,k));
for i:=k1 to n;
a2:=abs(a(k,i));
if a2<=a3 then l5;
i2:=i;
a3:=a2;
label l5;
next i;
if a1<=e and a3<=e then l2; for j:=k to n; label l6; for i:=1 to n-1; for i:=1 to n; label l2;
if i1:=k and i2=k then l3;
if a1
w:=a(i1,j);
a(i1,j)=a(k,j);
a(k,j):=w;
next j;
w:=b(i1);
b(i1):=b(k);
b(k):=w;
goto l3;
w:=a(k,k);
for j:=k to n;
a(k,j):=a(k,j)/w;
next j;
b(k):=b(k)/w;
for i:=1 to n;
if i=k then l7;
w:=a(i,k);
for j:=k to n;
a(i,j):=a(i,j)-w*a(k,j);
next j;
b(i):=b(i)-w*b(k);
label l7;
next i;
next k;
if x(i)=i then l8;
for j:=i+1 to n;
if x(j)<>i then l9;
w:=x(i);
x(i):=x(j);
x(j):=w;
w:=b(i);
b(i):=b(j);
b(j):=w;
j:=n;
label l9;
next j;
label l8;
next i;
print b(i);
next i;
end;
print "error";
end;
def fnf(x):=x*x-exp(-x);
input l,e;
label l1;
input a,b;
l1:=0;
if fnf(a)*fnf(b)>=0 l1;
label l2;
c:=(a+b)/2;
l1:=l1+1;
if fnf(a)*fnf(c)<0 then b:=c else a:=c;
if abs(fnf(c))>e and l1 print c,fnf(c); example:fnf(x):=x*x-exp(-x)=0
print l1;
end;
f)x):=0 > x:=f(x) >def
def fnf(x):=x^2/10+1;
input l,e;
input a;
c:=fnf(a);
for i:=1 to l;
c1:=c;
c:=fnf)c);
if abs(c-c1)
i:=i-1;
label l1;
print c,c-c1,i;
end;
def fnf(x):=3*-sin((x+y)/2);
fef fng(x):=3*y-cos((x-y)/2);
def fnh(x):=3-0.5*cos((x+y)/2);
fni(x):=-0.5*cos((x+y)/2);
fnj(x):=0.5*sin((x-y)/3);
fnk(x):=3-0.5*sin((x-y)/2);
input l,e;
input x,y;
for p:=1 to l;
f:=fnf(x);
g:=fng(x);
h:=fnh(x);
i:=fni(x);
j:=fnj(x);
k:=fnk(x);
if h*k-j*i=0 then l1;
x1:=(-f*k+g+i)/(h*k-j*i);
y1:=(-g*h+f*j)/(h*k-j*i);
x2:=x+x1;
y2:=y+y1;
if abs(x-x2)
y:=y2;
next p;
p:=p-1;
label l2;
print x2,x-x2,y2,y-y2;
print p;
label l1;
print "no converge";
end;
input n;
dim a(n+1),b(n+4),c(b+4);
input e;
for i:=n+1 to 1 step -1;
input a(i);
next i;
if n<3 then l1;
input p,q:
b(1):=0;
b(2):=1;
c(1):=0;
c(2):=1;
n:=n+1;
label l2;
for i:=1 to n-1;
k:=i+2;
b(k):=a(i+1)-p*b(k-1)-q*b(k-2);
c(k):=b(k)-p*c(k-1)-q*c(k-2);
print b(k),c(k);
next i;
n:=n+1;
c1:=c(n-1)-b(n-1);
d:=c(n-2)*c(n-2)-c1*c(n-3);
d1:=b(n-1)*c(n-2)-b(n)*c(n-3);
d2:=-b(n-1)*c1+b(n)*c(n-2);
d3:=d1/d;
d4:=d2/d;
n:=n-1;
print d,d1,d2,d3,d4;
p:=p+d3;
q:=q+d4;
if abs(d3/p)>e or abs(d4/q)>e then l2;
gosub l3;
n:=n-2;
if n:=2 then l4;
if n:=3 then l5;
for i:=2 to n;
a(i):=b(i+1);
next i;
goto l2;
label l5;
p:=b(3);
q:=b(4);
gosub l3;
end;
label l4;
print -b(4)/b(3);
end;
label l1;
if n=1 then print -a(1);:end;
p:=a(2);
q:=a(3);
gosub l3;
end;
label l3;
w:=p^2-4*q;
w1:=abs(w);
r:=sqr(w1);
if w<0 then l6;
r1:=0.5*(-p+r);
r2:=0.5*(-p-r);
i1:=0;
i2:=0;
goto l7;
label l6;
r1:=-0.5*p;
r2:=r1;
i1:=0.5*r;
i2:=-0.5*r;
label l7;
print r1,i1;
print r2,i2;
return;
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;