数値計算ツール・アラカルト

出典: 志村正人　シンポジウム　2004
論文:そるぶ(1）
PDF
Script
sorubu_806.ijs

関数の定義とグラフ

Example \( y=x^{3}+x-1 \)
Jでのp.を用いた定義。p. での定義では高次の項が右側に来る。　p.は小文字 
f01=: _1 1 0 1 &p.
サイエンスプロット: 書式 は plot　区間；'関数' 
 plot _5 5 ; 'f01'

非線形方程式の解法

Example: \( x^{4}-10x^{3}+36x^{2}-58x+35=0 \)
J の強力なpolynomial(p.)機能で解く 
   p. 35 _58 36 _10 1
+-+------------------------+
|1|4.41421 2j1 2j_1 1.58579|
+-+------------------------+

左の1は内部反復回数で苦労の目安。多いと苦労している

Newton法(一次）

Newton法/Newton-Raphson methodはJの簡潔なイディオムがある
　\[ x-\dfrac{f(x)}{f^{'}(x)} \]
Script New-11 は収束判定を行っているので希に失敗することがある。new-1は有限回ループ 
new_1=: 1 : ' ] - x % x D.1'(^:100) ("0)
new_11=: 1 : ' ] - x % x D.1'(^:_) ("0)
Example: \( f=x^{3}+x-1 \)
図で解の大まかな位置の見当をつける 
f=. _1 1 0 1&p.
   plot _3 3 ; 'f'
-3<-->3を初期値に選ぶ。Jはゾーンで選べる。 
   _1 1 0 1&p. new_1 i:3
0.682328 0.682328 0.682328 0.682328 0.682328 0.682328 0.682328

   p. _1 1 0 1
+-+---------------------------------------------+
|1|_0.341164j1.16154 _0.341164j_1.16154 0.682328|
+-+---------------------------------------------+

線形連立一次方程式

クラメール法

ガウスの吐き出し法が堅牢で定番だが、ここでは線形計算向きの優雅なクラメール法で解いてみよう。
Example: 連立方程式
\[ \left( \begin{array}{cccc} x_{1}& +2_{x2}& -x_{3}& = 2\\ &3_{x2}& +4_{x3}& = 18\\ x_{1}&& -x_{3}& = -2\\ \end{array} \right. \]
拡大係数行列 
   S1=: 3 4 $ 1 2 _1  2 0 3 4  18 1 0 _1 _2
   S1
1 2 _1  2
0 3  4 18
1 0 _1 _2
アルゴリズム
\[ \dfrac{拡大係数行列}{係数行列} \]
クラメイル法のスクリプト 
   cr=: %. }:"1
結果： 
   cr S1
1 0 0 1
0 1 0 2
0 0 1 3
  • 解は右の列
  • 左側が単位行列になればクリーンに解けている。
  • 数値計算のゴミが入ることがある。(numeric.ijsにあるcleanを先頭につけると取れる)
  • 左に他の数値が残るようなら行列が潰れており、うまく解けていない。

ガウス・ジョルダン法

addons/math/misc/linear.ijs

にJの線形計算ツールが入っている。 

   gauss_jordan   S1
1 0 0 1
0 1 0 2
0 0 1 3

非線形同時方程式

Example
\[ \left\{ \begin{array}{cc} f_{0}=&e^{x}+xy-1 \\ g_{0}=&sin xy + x + y +2\\ \end{array} \right. \]
関数の定義 
f00=: 3 : ’(ˆ{.y.)+(*/y.)-1’
g00=: 3 : ’2+(+/ y.)+1&o.*/ y.’
他変数のニュートン法の定義　17回まで反復 
new_2=: 1 : ’ ] - x. (%.|:)x. D.1’ (ˆ:17) ("1)
実行
　
   (f0,g0) new_2 1 1
_9.41145e_6 _2