load 'plot numeric trig'
load 'addons\lapack\lapack.ijs'
load 'addons\lapack\dgeev.ijs'

NB. -----------------------
cr=: %.}:"1 NB. Cramer     
NB. Newton
new_1=: 1 : ' ] - x. % x. D.1'(^:100) ("0)
new_11=: 1 : ' ] - x. % x. D.1'(^:_) ("0)
NB. Kerner
k=: 1 : '] - x. &p. % (<0 1)&|:@((1&(*/\."1))@(-/~))'
K=: k(^:_)
pp=:+//.@(*/)
cfr=:pp/@(- ,. 1:)
b=: 1 2 3 4
c=: cfr b+0.5
NB. diag
taikaku=: 3 : 0
NB. solve with disg
TMP0=: dgeev_jlapack_ y.
TMP1=: (%.>{: TMP0) +/ . * y. +/ . * >{:TMP0
)

new_2=: 1 : ' ] - x. (%.|:)x. D.1'  (^:17) ("1)


taikaku2=: 4 : 0
NB. x. is y_n
NB. seimal regression
TMP0=: dgeev_jlapack_ y.
TMP1=: (%.>{: TMP0) +/ . * y. +/ . * >{:TMP0
TMP2=: ( |: >{: TMP0) +/ . * x.
TMP3=: TMP2 %. TMP1
(> {: TMP0) +/ . * TMP3
)

  
NB. =======================
NB. Partial differential
NB. by C.Reiter
NB. ==========================
NB.name of  del f(x1 x2) is fixed  
a=: 3 5 2 *"1 (i.2 3 3 )e. 5 7 12
N=: 1{ $ a   NB. matrix size using in P
NB. P=: 1 : '((m. "_ +/ . * {:) +/ . * {.)@(^/&(i.3))'
P=: 1 : '((m. "_ +/ . * {:) +/ . * {.)@(^/&(i. N ))'

f=: a P

NB. jac=: (D.1)"1

NB. ===========================================
NB. Script is written by M. Shimura 08 July 2002
NB. ============================================  
new_ml=:  4 : 0
NB. Newton method for Maximum Likelyhood function
NB. x. Partial derivative matrix
NB. y. Start value
NB. N=: Matrix size 
Y=: y.
NR=: 0
while. NR < 10 do. NB. repeat 10 times (change if you wish)
D=: (%. x. P  D."1 Y) +/ . * x. P Y
Y=:NT=: Y - D
NR=: >: NR
end.
Y
)

NB. ====================
euler=: 1 : 0
NB. Usage: f euler y.
NB. u. is function //e.g. fn=: +/  is y'=y+x
NB. Usage: u. euler y. //e.g. f euler 1;0 1//(y0;y1)
NB. y0 is syokiti // y1 is  range of x  //e.g. 0 1 
BAND=:+/ | tmp=: > 1 { y.
dh=: BAND % 1024  NB.  dh //pitch is 1024
XSTEP=: ({. >{: y.) + (i.1024)  * dh
ANS=: ( (# XSTEP) , 2) $ 0
 NB. xn=: {.>{: y. NB. X0
 NB.  yn=: >{. y.
X=: {.>{: y. NB. X valuer
Y=: >{. y. NB. y1
COUNTER=: 0
while. COUNTER < # XSTEP do.
TMP=:dh *(TMP=: X,Y)
Y=:Y+ u. TMP NB. u. TMP is f(x,y)
ANS=:(X, Y) (COUNTER)}ANS
X=: X+dh
COUNTER=: >: COUNTER
end.
ANS
)


 runge=:1 : 0
 NB.  Runge Kutter Method
 NB.  y. A;B C //e.g. y.=: 2; 0 1
 NB. A is syokiti B C is start and goal of point to calc  
 NB. ussage: u. runge y. // u. is verb(function) ex. f=.+/ + */ 
 NB. e.g. fn=: (^&2@{.)%}. 
 BAND=:+/ | tmp=: > 1 { y.
 dh=: BAND % 1024  NB.  dh NB. pitch is 1024
 XSTEP=: ({. >{: y.) + (i.1024)  * dh
 ANS=: ( (# XSTEP) , 2) $ 0
  xn=: {.>{: y. NB. X0
  yn=: >{. y.
 COUNTER=:0
  while. COUNTER   < # XSTEP do.
NB.   if. COUNTER = PIECE do. goto_end. end. 
  x1=:xn + dh % 2
  x3=:xn + dh
  k0=: (u.  x3 , yn) * dh
  k1=: (u.  x1,y1=: yn +  k0 * 0.5) * dh
  k2=: (u.  x1,y2=: yn +  k1 * 0.5)* dh
  k3=: (u.  x3,y3=: yn + k2)* dh
  yy=:   (+/ k0 , (2 * k1 ) ,(2 * k2) , k3)  % 6
  yn=:yn + yy
  ANS=: ((COUNTER{XSTEP),yn) COUNTER } ANS
NB. ANS=:(  x3 , yn) ((COUNTER , 0) ;( COUNTER , 1)) } ANS  
  xn=:x3
  COUNTER=:>:COUNTER
  end.
NB. label_end.
  ANS
  )



NB. find Spectol
NB. ussage (\:m) sp matrix
m=:=/~@i.@#@[
n=:i.@#@[
tm=:|./(" 0 1) ~&n&[  NB. make index 
li=: [ */ m@[  NB. Lamda I
bb1=:([ tm ]) { [
bb2=: ( {.("1) - }.("1))@bb1
bb3=: (*/ "1) @ bb2
k1=:] - ("2) li NB. A-Lamda I
k2=: (<"2)@k1
k3=: >@ind { "1 k2
k4=:(+/ . * )/("3)
k5=:k4@(>@k3)
ind=: ( }. "1)@tm NB. index
sp=:k5 % bb3   NB. find spectol


spe=: 4 : 0
NB. spectol explicit
NB. x. is /: eigenvalue
NB. y. is matrix
NB. Usage. x. spe y. 
ind=: -. (i. # x.) =/ ~i. # x.
TMP0=: x. ,. ind #/ x.
TMP2=: */"1 TMP1=:({."1 TMP0) - }."1 TMP0 NB. bunshi
I=: =/~i. # y. NB. tani matrix
TMP3=: y. - L:0 <"2(}."1 TMP0) */ I
({@>% TMP2) * L:0  ({."1 TMP3) +/ . * L:0 {:"1 TMP3
if. 1=({: $ TMP3) do. ({@> % TMP2) * L:0 TMP3 end.
NB. if. double eigenbalue
)
 

NB. -------------------------------------------------
NB. multiple beki
NB. -------------------------------------------------
m_beki=: 4 : 0
NB. calc A^m
NB. x. is  m of A^m
NB. y. is matrix
TMP_0=: dgeev_jlapack_ y.
TMP_1=: (/:~>1{TMP_0) spe y. NB. find spectol
TMP_2=: ({@>/:~>1{TMP_0),. TMP_1
+/ > TMP_3=:(({."1 TMP_2)^ L:0  x.)* L:0 ({:"1 TMP_2)
NB.  usage: 1 m_beki KAS0
NB. 0 1
NB. 5 4
)


NB. -------define function--------------
f0=: 3 : ' y. - 2 o. y.'
f1=:_1 1 0 1 &p.

f0=: -&1@(*/)+^@{.
g0=: +&2@(+/)+1&o.@*/

f00=: 3 : '(^{.y.)+(*/y.)-1'
g00=: 3 :  '2+(+/ y.)+1&o.*/ y.'


jac0=: 1 : 'x. D.1'

NB. ------------data block----------------
S0=: 3 4 $ 1 4 3 1 2 5 4 4 1 _3 _2 5
S1=: 3 4 $ 1 2 _1  2 0 3 4  18 1 0 _1 _2
S2=: 3 3 $ 0 1 0 0 0 1 _0.25 0.25 1
S3=: 3 3 $  6 4 _2 4 12 _4 _2 _4 13
S12=: 1 3 3 1
S13=: 0 1 2 3 4 5

a7=:3 3 $ 1 1 _1 1 1 _1 _1 1 1
 
a8=: 3 3 $ 0 1 1 1 0 1 1 1 0
a9=: 3 3 $ 1 2 1 _1 4 1 2 _4 0