"rootSolve"の中の multiroot() 多元の連立方程式の解（のベクトル）を返してくれる library(rootSolve) model <- function(x) c(F1= x[1] + x[2] + x[3]^2 - 12, F2= x[1]^2 - x[2] + x[3] - 2, F3= 2 * x[1] - x[2]^2 + x[3] - 1 ) (ss<-multiroot(model,c…

こちらで連立方程式の解を求めたい パッケージ"rootSolve" の中の uniroot.all() で解を出す この uniroot.all() はパッケージ"base" の中の uniroot() 関数を使っている uniroot()の中身の一部 val <- .Internal(zeroin2(function(arg) f(arg, ...), lower,…

線型連立方程式を解く 以前はこちらやこちら 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,]…

こちらとこちらの続き 線維方向を考えたうえで3Dに拡張する 6近傍への移動を考える 3Dで線維方向を含めた誘電率の値はこちらに 頂点への移動が少ない可能性がある library(rgl) Nx<-30 Ny<-30 Nz<-30 Nt<-100 U<-array(0,c(Nx,Ny,Nz)) U[11:20,11:20,11:20]<…

こちらのつづき kの値を各位置について与える まず を各位置で与え、その から k を計算する Nx<-30 Ny<-30 Nt<-100 U<-tempU<-matrix(0,Nx,Ny) U[11:20,11:20]<-1 #初期分布 #これがl方向とt方向の拡散の定数 kl<-0.20 kt<-0.05 #線維走向の傾き theta<-mat…

非線形連立方程式を解きたい こちらの中のこちらのPDFを参考にして まず二分法で を解いてみる f <- function(x){ sin(x)-1/2 } d<-1e-15 xl<-0 xr<-pi/2 xm<-(xl+xr)/2 while(abs(f(xm))>d){ if(f(xl)*f(xm)>0){ xl<-xm }else if(f(xl)*f(xm)<0){ xr<-xm } …

ナナメ方向への拡散をあらわす こちらのつづき 角度を変えて試してみると、斜めの拡散は縦横の移動分とそれに対する45度方向(kxy,kyx)の移動が合わさって実現されてるようだ 縦横の成分 -dUx*kxx, -dUy*kyy を 0 にして斜めにだけ移動させて計算すると、いつ…