スタートページJavaScriptライブラリ目次

数値解析関連 JavaScript関数ライブラリ na.js の利用解説書

ご利用にあたって

代数系  f(x) を指定する


代数方程式 f(x) = 0 の根:はさみうち法
  naFx0HasamiuchiTrace(表示場所, fx, xmin, xmax[, eps]) 
  var x = naFx0Hasamiuchi(fx, xmin, xmax[, eps]) 

はさみうち法による f(x) = 0 の解(参照:「方程式の根」)を求める。
fxは文字列として、全体を" "で囲む。
xmin~xmaxの間に唯一の実根があることを前提にしている( f(xmin) と f(xmax) が同符号ではいけません。
二分法でxmin~xmaxの範囲を狭めていき、絶対誤差 |f(x)| < eps になるまで繰り返す。epsを省略すると、(xmax-xmin)/1000とする。

ケース1

入力
  var fx = "x*x-1";
  var xmin = 0;
  var xmax = 3;
  var eps = 0.001;
出力
  var x = naFx0Hasamiuchi(fx, xmin, xmax, eps);
//  x:1

ケース2

fx= xmin= xmax= eps=

naFx0HasamiuchiPl 代数方程式 f(x) = 0 の複数の実根:はさみうち法
  naFx0HasamiuchiPlTrace(表示場所, fx, xmin, xmax) 
  var x = naFx0HasamiuchiPl(fx, xmin, xmax) 

はさみうち法による f(x) = 0 の1個の実数解を求める関数にはnaFx0Hasamiuchiがあります。
それに対して、naFx0HasamiuchiPlTracePl は、複数の実数解を求める関数です。重根は1個の根とします。
xmin~xmax の範囲を大雑把にスキャンして、 f(x) の正負が逆転する範囲 x~x+dx の組を探索し、dmin(k) = x, dmax = x+dx とします。
そして、それぞれの dmin~dmax について外部関数 naFx0Hasamiuchi を呼び出しています。
コードをダウンロードするときには、naFx0Hasamiuchi もダウンロードしてください。
本当は根があるのに見落としたり、NaN になる場合があります。

ケース1

入力
  var fx = "x*(x-1)*(x+1)";
  var xmin = -2;
  var xmax = 2;
  var x = naFx0HasamiuchiPl(fx, xmin, xmax);
出力
  x:[-1, 0, 1]

ケース2

fx= xmin= xmax=

方程式の根(ニュートン法)
  naFx0NewtonTrace(表示場所, fx, x0[, eps]);  
  var x = naFx0Newton(fx, x0[, eps]); 

ニュートン法により、方程式の実根を求める。
   xn = xn-1 - f(xn-1)/f'(xn-1)   (f'=df/dx)
を繰り返して根に近づく。参照:「方程式の根」
 x0は初期値。その値により異なる解に達したり、解が得られないことがある。

fxは文字列として、" "で囲む。eps を省略すると 0.001 になる、

ケース1

入力
  var fx = "x*x-1";
  var x0 = 0.3;
  var eps = 0.001;
出力
  var x = naFx0Newton(fx, x0, eps);
//  x:1

ケース2

fx= x0= eps=

naFx0NewtonPl 代数方程式 f(x) = 0 の複数の実根:ニュートン法
  naFx0NewtonPlTrace(表示場所, fx, xmin, xmax) 
  var x = naFx0NewtonPl(fx, xmin, xmax) 

ニュートン法による f(x) = 0 の1個の実数解を求める関数にはnaFx0Newtonがあります。
それに対して、naFx0NewtonPlTracePl は、複数の実数解を求める関数です。重根は1個の根とします。
xmin~xmax の範囲を大雑把にスキャンして、 f(x) の正負が逆転する範囲 x~x+dx の組を探索し、x0(k) = x とします。
そして、それぞれの k について、初期値を x0(k) として外部関数 naFx0Newton を呼び出し、解 x[k] を求めます。
コードをダウンロードするときには、naFx0Newton もダウンロードしてください。
本当は根があるのに見落としたり、NaN になる場合があります。

ケース1

入力
  var fx = "x*(x-1)*(x+1)";
  var xmin = -2;
  var xmax = 2;
  var x = naFx0NewtonPl(fx, xmin, xmax);
出力
  x:[-1, 0, 1]

ケース2

fx= xmin= xmax=

3次方程式の複素数根 a0 + a1x + a2x2 + a3x3 = 0
  naFx03jiTrace(表示場所, a0, a1[, a2[, a3]])  
  var x = naFx03ji(a0, a1[, a2[, a3]])   

 

a3を省略すると2次方程式、さらにa2を省略すると1次方程式になる。
2次方程式は根の公式による。3次方程式は、特殊なケースを除き、ニュートン法により一つの実根を求め、式の変形により2次方程式にして求める。
戻り値は2次元配列である。下の[出力]参照

ケース1

入力
  var x = naFx03ji(2, 0, -1, 1); // 2x3 - x + 1 = 0
出力
  x[k][0] 実数部 x[k][1] 虚数部
  x[0][0]= 1.000 x[0][1]= 0.9999
  x[1][0]= 1.000 x[1][1]=-0.9999
  x[2][0]=-1.000 x[2][1]= 0

ケース2

a0 + a1x + a2x2 + a3x3 = 0

数式 f(x) を最大にする x
  naFxMaxTrace(表示場所, fx, xmin, xmax [, eps])  
  var x = naFxMax(fx, xmin, xmax [, eps]) 

数式 f(x) が最大となるxの値を戻します。
 xmin~xmaxを100等分してdxとし、for (x=xmin, x<=xmax; x=x+dx) でf(x) を最大にする x を求め、次に x-dx~x+dx で同様な操作を繰り返すことにより、最終的な最適値 x を求めます。epsを省略すると(xmax-xmin)/1000になります。

ケース1

入力
  var fx= "x*(x-1)*(x-2)";
  var xmin = -1;
  var xmax = 2;
  var eps = 0.0001;
出力
  var x = naFxMax(fx, xmin, xmax, eps);
  // x = 0.4226

ケース2

fx= xmin= xmax= eps=

代数式 f(x) のxにおける微分係数a
  naFxDiffTrase(表示場所, fx, x [, d, eps])  
  var a = naFxDiff(fx, x [, d, eps])  

fxのxにおける微係数aを求めます。xの近傍の2点 x-d*xとx-d*xでの勾配をa0とし、dを1/2にしたときの勾配aとの差が、相対誤差 |(a-a0)/a| < eps になるまで繰り返します。
dを省略すると0.01、epsを省略すると0.001にしています。

ケース1

入力
  var fx= "x*x";
  var x = 1;
  var d = 0.01;    // 微小区間
  var eps = 0.001;  // 打切誤差
  var a = naFxDiff(fx, x, d, eps);
出力
  a = 2

ケース2

fx= x= d= eps=

微分方程式(ルンゲクッタ法)
  naFxRungeKutta(表示場所, fxy, x, dx)  
  a = naFxRungeKutta(fxy, x, dx)  

一階微分 dy/dx = f(x,y) の式を与えて y(x+dx) = y(x) + a*dx の a の値を、4次のルンゲクッタ法により求めます。

ルンゲクッタ法では、次のステップで定義します。
    k1 = f(x,y)
    k2 = f(x+dx/2, y+k1*dx/2)
    k3 = f(x+dx/2, y+k2*dx/2)
    k4 = f(x+dx, y+k3*dx)
   a = (k1 + 2*k2 + 2*k3 + k4) / 6

ケース1

入力
  var fxy= "x*y";  // 文字列 f(x,y) をJavaScript形式で記述
  var x = 2;    // fxy に x がなくても形式的に適当な値を与える
  var y = 1;    // fxy に y がなくても形式的に適当な値を与える
  var dx = 0.01;
  var a = naFxRungeKutta(yx, x, y, dx);
出力
  a=2.0252

ケース2

fxy= x= y= dx=

微分方程式(ルンゲクッタ法)作表
  var [x, y, a, dy] = naFxRungeKuttaTable(fxy, xmin, xmax, dx, y0)  

一階微分 dy/dx = f(x,y) の式を与えて 4次のルンゲクッタ法により x = xmin~xmax に対応する y の値、および関連する値を算出します(y0 は x = xmin のときの y の値)。
  x: x = xmin, xmin+dx, xmin+2*dx, …, xmax の値
  y: x に対応する y の値
  a: (x,y) における微分係数
  dy: y の増分 a*dx
ここでは参考のために、表の一部を表示していますが、関数自体は表示機能は持っていません。

dy/dx=fxy= xmin= xmax= dx= y0=

連立微分方程式(ルンゲクッタ法)
  naFxRungeKutta2Trace(表示場所, ft, gt, x, y, dt)  
  var [ax, ay] = naFxRungeKutta2(ft, gt, x, y, dt)  

二つの一階微分 dx/dt = f(x,y), dy/dt = g(x,y) の式を与えて、
点(x,y) で x(t+dt) = x(t) + ax*dt, y(t+dt) = y(t) + ay*dt となる偏微分係数 ax, ay の値を
4次のルンゲクッタ法により求めます。

ルンゲクッタ法では、次のステップで定義します。
    k1 = f(x,y)
    k2 = f(x+dx/2, y+k1*dx/2)
    k3 = f(x+dx/2, y+k2*dx/2)
    k4 = f(x+dx, y+k3*dx)
   a = (k1 + 2*k2 + 2*k3 + k4) / 6

ケース1

入力
  var ft = "4*x - x*y";    // dx/dt を 文字列 ft として与える
  var gt = "-10*x + 2*x*y"; // dy/dt を 文字列 gt として与える
  var x = 10;
  var y = 15;
  var dt = 0.01;
  var ax. ay = naFxRungeKutta2(ft, gt, x, y, dt);
出力   ax = -104.728, ay = 221.522

ケース2

dx/dt=ft= dy/dt=gt= x= y= dt=

連立微分方程式(ルンゲクッタ法)作表
  var [t, x, ax, dx, y, ay, dy] = naFxRungeKutta2Table(ft, gt, x0, y0, tmin, tmax, dt);

二つの一階微分 dx/dt = f(x,y), dy/dt = g(x,y) の式と t = tmin における x0, y0 を与えて、t = tmin~tmax (刻み dt)での x(t), y(t) などを表の形式で戻します。

  t: t[i] = t t = tmin, tmin+dt, tmin+2*dt, …, tmax の値
  x: x[i]
  ax:偏微分係数 ∂x/∂t(t)
  dx:x の増分 ax[i]x*dt
  y: y[i]値
  ay:偏微分係数 ∂y/∂t(t)
  dy:y の増分 ay[i]*dx
ここでは参考のために、表の一部を表示していますが、関数自体は表示機能は持っていません。

dx/dt=ft= dy/dt=gt= x0= y0= tmin= tmax= dt=

数式 f(x) の極値(複数)
  naFxExtremeTrace(表示場所, fx, xmin, xmax [, eps])  
  var [極値, 大小] = naFxExtreme(fx, xmin, xmax [, eps]) 

数式 f(x) の極値 f'(x) が正から負へ、負から正に変化する点(複数)を戻します。
 まず、xmin~xmaxを100等分してdxとし、x-dx と x+dx の間で f'(x) の符号が変化したら極値があることになります。
 次に x-dx~x+dx で同様な操作を繰り返すことにより、最終的な最適値 x を求めます。精度は0.0001 になります。
戻り値の形式
 極値:xi
 大小:xi が極大なら1、極小なら-1

ケース1

入力   var fx= "x*(x-1)*(x-2)";
  var xmin = -1;
  var xmax = 3;
  var [極値, 大小] = naFxExtreme(fx, xmin, xmax);
出力
  極値[0] = 0.423 大小[0] = 1  // x=0.423 のとき、極大になる
  極値[1] = 1.578 大小[1] = -1 // x=1.578 のとき、極小になる

ケース2

fx= xmin= xmax=

代数式 f(x) の定積分(シンプソン則)
  naFxSimpsonTrace(表示場所, fx, xmin, xmax [, n2, eps]) 
  var s = naFxSimpson(fx, xmin, xmax [, n2, eps]) 

シンプソン則は、2次式で近似します。
 閉区間[xmin,xmax]を2n等分して、h=(xmax-xmin)/2にすると(すなわち
   x0(=xmin),x1,x2,....,x2n(=max)
ととして、xiに対応するfxの値を yiとすると)
   S=(xmax-xmin)/6n *{(y0+y2n)+4(y1+y3+...y2n-1)+2(y2+y4+...y2n-2) }
となります。参照:「定積分」

最初のn2で計算した値をs0、次にn2を2倍(すなわちキザミ幅を1/2)にして計算した値を s として、相対誤差 |(s-s0)/s| <eps となるまで繰り返します。
n2を省略すると20、epsを省略すると (xmax-xmin)/1000 にしています。

ケース1

入力
  var fx = "x*x";
  var xmin = 1;
  var xmax = 2;
  var n2 = 20;
  var eps = 0.001;
出力
  var s = naFxSimpson(fx, xmin, xmax, n2, eps);
  // s = 2.333

ケース2

fx= xmin= xmax= n2(キザミ数)= eps=

漸化式の特性方程式

漸化式の係数と初期値を与えて、an の一般解を「文字列」として戻します。
ここでは一般解が1命令で表現できる漸化式に限定しています。
Trace では、元の漸化式を順次計算した数列と、一般解にnを代入した数列を表示するので、一般解の確認ができます。
(係数と初期値の値によっては、一般解が得られないこともあります)

 an+1=p*xn+q
  naZenka2Trace(表示場所, a1, p, q) 
  var an = naZenka2(a1, p, q) 

  a1=  p= q= an=

 an+1=p*xn+ q + rn
  naZenka2nTrace(表示場所, a1, p, q, r) 
  var an = naZenka2n(a1, p, q, r) 

  a1=  p= q= r= an=

 an+1=p*an+ q*rn
  naZenka2eTrace(表示場所, a1, p, q, r) 
  var an = naZenka2e(a1, p, q, r) 

  a1=  p= q= r= an=

 an+1=p*anq
  naZenka2expTrace(表示場所, a1, p, q) 
  var xn = naZenka2exp(a1, p, q) 

a1 および(あるいは)p が負数のときは、q を整数に限定しています。
その場合でも、適切な解がえられないことがあります、
  a1=  p= q= an=

 an+2=p*an+1+ q*an
  naZenka3Trace(表示場所, a1, a2, p, q) 
  var an = naZenka3(a0, a1, p, q) 

特性方程式の判別式 D = q2 + 4q の正負により、計算方法が異なります。
D > 0 (実根)の例:a1=0, a2=1, p=2, q=3
D = 0 (重根)の例:a1=0, a2=3, p=6. q=-9
D < 0 (虚根)の例:a1=2, a2=4, p=2, q=-3
  a1= a2=  p=  q= an=

多点系  点(x,y) を与える


多点系 最小二乗法


最小二乗法(1次式)y = a[0] + a[1]x
  naLsm1Trase(表示場所,x,y)
  var a = naLsm1(x,y)

ケース1

入力
  var x = [ 0, 2, 4, 6, 8, 10, 12 ];
  var y = [ 1, 20, 30, 45, 50, 60, 70 ];
  var a = naLsm1(x, y);
出力
  a[0] = 6.54, a[1] = 5.48
  // y = 6.54 + 5.48x

ケース2

x=
y=

最小二乗法(1次式)の拡張 相関係数と95%信頼区間
  naLsm1ExTrace(表示場所, x,y)
  var [a, r, CI] = naLsm1Ex(x,y)

ケース1

入力
  var x = [ 0, 2, 4, 6, 8, 10, 12 ];
  var y = [ 5, 20, 30, 45, 50, 60, 70 ];
  var rtn = naLsm1Ex(x, y);
出力
  係数 a[0]=8.393, a[1]=5.268
  相関係数 r=0.993
  信頼区間 CI[0]=4.830, CI[1]=3.789, … ,CI[6]=4.830

ケース2

x=
y=

最小二乗法(n次式)y = a[0] + a[1]x + … + a[n]*xn の a[i] を求める
  naLsmMDTrace(表示場所, x,y.n)
  var a = naLsmMD(x,y.n)

n次式、
   y=a[0]+a[1]x+a[2]x2+・・・+a[n]xn
の a[i] は、次の連立方程式の解になります。

a[0]a[1]a[2]・・・a[n]定数項
n∑x∑x2・・・∑xn∑y
∑x1∑x2∑x3・・・∑xn+1∑yx
∑x2∑x3∑x4・・・∑xn+2∑yx2
∑xn∑xn+1∑xn+2・・・∑x2n∑yxn

ケース1

入力
  var y = [ 8, 5, 2, 0, 3, 6, 8 ];
  var x = [-3, -2, -1, 0, 1, 2, 4 ];
  var n = 2;
  var a = naLsmMD(x, y, n);
出力
  a[0] = 2.12, a[1] = -0.27, a[2]=0.50
  // y = 2.12 - 0.27 x + 0.50 x2

ケース2

y=
x=
n=


最小二乗法(多項式)y = a[0] + a[1]x[0] + … + a[n]*x[n-1] の a[i] を求める
  naLsmMVTrace(表示場所, x,y)
  var a = naLsmMV(x,y)

多項一次式
   y=a[0]+a[1]x[0]+a[2]x[1]+・・・+a[n]x[n-1]
の係数 a[i] は、次の連立方程式の解になります。

 a[0] a[1] a[2] ・・・ a[n]  定数項
 n ∑x[0] ∑x[1] ・・・ ∑x[n-1]  ∑y
 ∑x[0] ∑x[0]x[0] ∑x[0]x[1] ・・・ ∑x[0]x[n-1]  ∑yx[0]
 ∑x[1] ∑x[1]x[0] ∑x[1]x[1] ・・・ ∑x[1]x[n-1]  ∑yx[1]
 : : : : :  :
 ∑x[n-1] ∑x[n-1]x[0] ∑x[n-1]x[1] ・・・ ∑x[n-1]x[n-1]  ∑yx[n-1]

ケース1

入力
  var y =  [12, 15, 4, 2, 18, 1, 14];
  var x = [ [ 0, 1, 2, 3, 4, 5, 6],  // x[0][j]
       [ 2, 5, 3, 1, 7, 1, 4],  // x[1][j]
       [ 4, 3, 0, 1, 5, 2, 6] ];  // x[2][j]
  var a = naLsmMVTrace(x, y);
出力
  a[0]=1.08, a[1]=-1.046, a[2]=1.68, a[3]=1.98
  // y= 1.08 - 1.04 x[0]+1.68 x[1]+1.98 x[2]

ケース2

  y=
x[0]=
x[1]=
x[2]=
x[3]=


多点系 補間法


ニュートンの補間法
  naPolNewtonTrace(表示場所, x, y, px) 
  var py = naPolNewton(x, y, px) 

n+1個の点(x0,y0), (x1,y1), …,(xn,yn)を通るn次式 f(x)を、
   f(x)=c0 + c1(x-x0) + c2(x-x0)(x-x1) + … + cn(x-x0)(x-x1)…(x-xn-1
の形で表します。そして、
   f(x0)=y0、f(x1)=y1、…、f(xn+1)=yn+1
の関係から、cの値を求めて、n次式 f(x) を決定します。参照:「ニュートンの補間公式」

その式に、x=px に対応する y の値 py を戻します。

ケース1

入力
  var x = [ 0, 2, 3, 5 ];
  var y = [ 1, 5, 8,10 ];
  var px = [1, 4];
  var py = naPolNewton(x, y, px);
出力
  py[0] = 2.267, py[1]=10.067

ケース2

x[i]
y[i]
px[k]

ラグランジュの補間法
  naPolLagrangeTrace(表示場所, x, y, px) 
  var py = naPolLagrange(x, y, px) 

n+1組の点xa[0]~xa[n]、ya[0]~ya[n] を通るn次式をラグランジュの補間法により求め、m 個の、xa[n+1]~xa[n+1} における ya[n+1]~ya[n+m] を戻します。

n次式に一般化して、
   y=c0y0+c1y1+c2y2+・・・+cnyn
としたときの、ckは、次の式になります。
  ck = (x-x0)/(xk-x0) * (x-x1)/(xk-x1) * … * (x-xk-1)/(xk-xk-1)
          ((x-xk)/(xk-xk)の項がないことに注意)
    * (x-xk+1)/(xk-xk+1) * … * (x-xn)/(xk-xn)

 この式に、x = px を与えたときの y の値 py を戻します。
参照:「ラグランジュの補間法」

ケース1

入力
  var x = [ 0, 2, 3, 5, 7, 8 ];
  var y = [ 4, 3, 6, 8, 6, 8 ];
  var px = [1, 4. 6 ];
  py = naPolLagrange(x, y, px);
出力
  py = 1.067, 7.905, 6.829

ケース2

x[i]
y[i]
px[k]

スプライン補間(3次)
  naPolSplineTrace(表示場所, x, y, px) 
  var py = naPolSpline(x, y, px) 

スプライン補間は、n+1個の点 (x0,y0), (x1,y1), … , (xn,yn) を、2つの点の区間 [xj~xj+1]ごとに別々の式(一般には3次式)Sj(x) を設定して緩やかな曲線にし、極端な補間値にならないようにする補間法です。参照:「スプライン補間」

区間 [xj~xj+1] における3次式を
  Sj=aj(x-xj)3 + bj(x-xj)2 + cj(x-xj) + dj   (j=0~n-1)
とします。
ここで aj, bj, cj, dj を求めることになります。

これらの3次式が滑らかにつながっている必要があります。それで、xj におけるSj-1とSjの1次導関数(傾き)S(1) と2次導関数(滑らかさ)S(2) の値が等しいという条件を加えます。さらに、始点 x0 と終点 xn での2次導関数の値が0であるとの仮定をします。

ケース1

入力
  var x = [ 0, 2, 3, 5, 7, 8 ];
  var y = [ 4, 2, 6, 8, 6, 8 ];
  var px = [1, 4. 6 ];
  var py = naPolSpline(x, y, px);
出力
  py = 1.614, 8.189, 6.499

ケース2

x[i]
y[i]
px[k]

最小二乗法(高次式)による補間
  naPolLeastSquareTrace(表示場所, 次数, x, y, px> 
  var [c, py] = naPolLeastSquare(次数, x, y, px>) 

n組の(x,y)から、最小二乗法により、r次の近似式 y=c0+a1x+c2x2+...+crxr の係数 ci を求め。この近似式により、f個の px に対応する py を計算します。;
参照:「回帰分析・最小二乗法」「多項式、高次式での最小二乗法」

ケース1

入力
  var 次数 = 3;
  var x = [ 0, 2, 3, 5, 7, 8 ];
  var y = [ 4, 2, 6, 8, 6, 8 ];
  var px = [1, 4. 6 ];
  var [c, py] = naPolLeastSquare(表示場所, 次数, x, y, px);
出力
  c = 3.610, -0.778 ,0.523, -0.046
   // 近似式 = 3.610 - 0.778x + 0.523 x2 - 0.046 x3
  py = 3.309, 5.894, 7.731

ケース2

次数
x[i]
y[i]
px[k]

ベジェ曲線(2次、3次)
  var [X, Y] = naBezier(x, y, p) 

ベジェ曲線とは、始点と終点を固定し、その直線を制御点から引っ張るような曲線にしたものです。
その計算方法は下式の通りです。
これらの式で、t=0 ~ 1 を p個に等分して、曲線の座標を (X[i], Y[i]) i=0 ~ p で戻します。
 2次ベジェ(3点を与える。始点(x[0,y[0]), 終点(x[2],y[2]), 制御点(x[1],y[1])
  X[i] = (1-t)2*x[0] + 2(1-t)t*x[1] + t2*x[2]
  Y[i] = (1-t)2*y[0] + 2(1-t)t*y[1] + t2*y[2]
 3次ベジェ(4点を与える。始点(x[0,y[0]), 終点(x[3],y[3]), 制御点(x[1],y[1]),(x[2],y[2])(x[2],y[2])
  X[i] = (1-t)3*x[0] + 3(1-t)2t*x[1] + 3(1-t)*t2*x[2] + t3*x[3]
  Y[i] = (1-t)3*y[0] + 3(1-t)2t*y[1] + 3(1-t)*t2*y[2] + t3*y[3]
式を求める手順は単純ですので、trace は省略します。

ケース1

      始点  制御点 終点
  var x = [ 0, 1,  2,  3];
  var y = [ 0, 0.5, 3.5, 3];
                p
  var [X, Y] = naBezier(x, y, 10);

ケース2

x[0]= y[0]=
x[1]= y[1]=
x[2]= y[2]=
x[3]= y[3]=
p==10,50,100

   

正方行列系 線形連立方程式など


連立一次方程式(ガウス=ジョルダンの掃出し法)
  naGaussJordanTrace(表示場所, 係数表, 定数項) 
  var x = naGaussJordan(係数表, 定数項) 

説明変数や被説明変数の標本から、連立方程式の定式化を行い、解を求めるまでの一連の処理を行うには、 maMregEq関数があります。

ケース1

入力
  var 係数表 = [ [0, 3, 2],    // 0x0 + 3x1 + 2x2 = 12
          [2, 4, 0],    // 2x0 + 4x1 + 0x2 = 10
          [3, 1, 5] ];   // 3x0 + 1x1 + 5x2 = 20
  var 定数項 = [12, 10, 20];
  var x = naGaussJordan(係数表, 定数項);
出力
  x[0] = 1, x[1] = 2, x[2] = 3

ケース2

定数項 x0 x1 x2 x3 x4
0
1
2
3
4

naTdma 三重対角行列の連立方程式
  naTdmaTrace(表示場所, 三重対角行列);  
  var x = naTdma(三重対角行列);  

三重対角行列とは、下図の左行列のように、対角線とその左右要素以外の要素が0である行列です。右のd行列を定数項としたときの解xを求めます。
1次元の移流方程式や拡散方程式などを偏微分方程式は、このような三重対角行列になることが多くあります。
ガウスの消去法よりも効率的な計算方法にTDMA(TriDiagonal-Matrix Algorithm):三重対角行列アルゴリズム)があります。

三重対角行列はゼロ要素が多く、入力データを行列形式で与えるのは非効率ですので、ここでは次の形式で与えることにします。
  var a = [a0, a1, …, an-1]
  var b = [b0, b1, …, bn-1]
  var c = [c0, c1, …, cn-1]
  var d = [d0, d1, …, bn-1]
   (a0とcn-1は不要なのですが、形式を整えるために、任意の値(0)を与えることにします。)

ケース1

モデル
  三重対角行列(係数)  定数項
  ┌ 1 3       ┐  ┌ 7 ┐
  │ 2 2 1     │  │ 7 │
  │   1 1 2   │  │ 3 │
  │     3 2 1 │  │ 5 │
  └       3 3 ┘  └ 6 ┘
入力データ
  var a = [0, 2, 1. 3, 3];
  var b = [1, 2, 1, 2, 3];
  var c = [3, 1, 2, 1, 0];
  var d = [7, 7, 3, 5, 6];
出力
  x[0] = 1;
  x[1] = 2;
  x[2] = 1;
  x[3] = 0;
  x[4] = 2;
検算
  0行: 1*1 + 3*2          = 7
  1行: 2*1 + 2*2 + 1*1       = 7
  2行:    1*2 + 1*1 + 2*0    = 3
  3行:       3*1 + 2*0 + 1*2 = 5
  4行:          3*0 + 3*2 = 6

ケース2

 a b c d

N+1個の点を通るN次式の係数
  naCoefNthDegTrace(表示場所, x, y) 
  var c = naCoefNthDeg(x, y) 

例えば、4(=N+1)個の点 [x0,y0], [x1,y1], [x2,y2], [x3,y3] を通る3(=N)次式
   y = c[0] + c[1] x + c[2] x2 + c[3] x3
の係数 c[i] を求めます。
x0 < x1 < … < xn にしてください。xn が極端に大きくならないように単位を調整してください。
連立一次方程式(ガウス=ジョルダンの掃出し法)の応用です。

ケース1

入力
  var x = [-1, 0, 1, 2];  // 4つの点→3次式
  var y = [ 2, 2, 6, 8];
  var c = naCoefNthDeg(x, y);
出力
  c[0] = 2, c[1] = 3, c[2] = 2, c[3] = -1
  // y = 2 + 3 x + 2 x2 - 1 x3

ケース2

x=
y=