状態の変化を表した式から平衡状態を求めたい 1次元のブリュッセレータで試してみる パッケージrootSolveの中のmultiroot関数を使う multirootの中身でもRの逆行列（solve）をつかっているようだ 黒がX、赤がYの濃度 解が安定かどうかはたぶんまた別のはなし…

連立方程式の計算 一般化、大規模化のためにこちらの改良 Rの逆行列を用いている x<-c(0,0) T<-100 d<-0.01 dim<-length(x) f<-function(x){ z<-rep(0,dim) z[1]<-(x[1]-3)^2+x[2]^2-3 z[2]<-sin(x[1])+exp(x[2]-1)-1 return(z) } g<-function(i){ v<-rep(0,…

ローレンツ方程式（Wiki） 位相次元 フラクタル次元 フラクタル次元の中に、情報次元、ボックスカウント次元、ハウスドルフ次元、レニー次元などの定義があるらしい # ローレンツ方程式 p<-10 r<-28 b<-8/3 f<-function(x,y,z){-p*x+p*y} g<-function(x,y,z)…

非線形連立方程式 解の数値計算の方法 こちら のpdf より となるように を選ぶとすると、2次より高次の項を無視すると を得る を繰り返し計算する 以前の rootSolve の記事 f<-function(x,y){(x-3)^2+y^2-3} g<-function(x,y){sin(x)+exp(y-1)-1} # fx<-func…

Solving Equation on the Computer と こちら のコメントで教えていただいた数値計算について 数値計算 Euler's method Improved Euler method Runge-Kutta method その他（Wiki） について Euler's method Improved Euler method 2点の傾きの平均をとる Run…

線型連立方程式を解く 以前はこちらやこちら Gauss-Seidel 法 , は の対角成分からなる対角行列 対角成分が大きくないといけないらしい N<-3 A<-matrix(runif(N^2),N,N) diag(A)<-diag(A)*10 #対角成分が大きくないといけない x<-runif(N) b<-runif(N) D<-ma…

こちらのつづき ピボット操作というのは入っていないが、掃き出し法での式変形をもう少し再現できるようにしてみる N<-4 M<-5 A<-matrix(runif(N*M),N,M) C<-A L<-min(N,M-1) for(i in 1:L){ C[i,]<-C[i,]/C[i,i] for(j in 1:N){ if(j != i){ C[j,]<-C[j,]-C…

パッケージに "rootSolve"というものがあるらしい こちらのコメントで教えていただいたもの 詳しくは次を実行 vignette("rootSolve") 中身は以下の通り rootSolve uniroot.all : to solve for all roots of one (nonlinear) equation multiroot : to solve n…

線形の連立方程式の解法 こちらのつづき 掃き出し法に近いもので、Gauss-Jordan法 拡大係数行列を変形させてみる N<-3 A<-matrix(runif(N*(N+1)),N,N+1) C<-A for(i in 1:N){ C[i,]<-C[i,]/C[i,i] if(i != N){ for(j in (i+1):N){ C[j,]<-C[j,]-C[j,i]*C[i,]…