# 質問部屋
```julia
using Images
using FFTW
using NLsolve
#テストするための画像を作る．
u = zeros(128,128) #真っ黒な画像
# make a view
img_a = @view u[1:128,24:48] #この行列を変化させると大元の行列も変えることができる．
img_a = fill!(img_a, 1.0)
img_b = @view u[1:128,72:96]
img_b = fill!(img_b, 1.0)
img_c = @view u[1:128,104:112]
img_c = fill!(img_c, 1.0)
img_d = @view u[40:70,1:128]
img_d = fill!(img_d, 0.5)
f = u #動かしたくない領域を参照するための行列
u #動かすための行列．こいつを更新していく
dt = 0.05
e = 0.1
C1 = 300
C2 = 150000
lam = 50000
function nabla1(u)
fftu = fft(u)
fftu[1,:] = fftu[1,: ] * 0
for j in 1:63
fftu[j+1,: ] = (2 *pi *j *1im /128 ) * fftu[j+1 ,: ]
fftu[129-j ,: ] = (2 *pi *(-j) *1im /128 )*fftu[129-j ,: ]
end #j
real(ifft(fftu))
end
function nabla2(u)
fftu = fft(u)
fftu[: ,1 ] = fftu[: ,1 ] * 0
for k in 1:63
fftu[: ,k+1 ] = (2 *pi *k *1im /128 ) * fftu[: ,k+1 ]
fftu[: ,129-k ] = (2 *pi *(-k) *1im /128 )*fftu[: ,129-k ]
end #k
real(ifft(fftu))
end
function lap(u) #ラプラシアンの微分演算子．
nabla1(nabla1(u)) + nabla2(nabla2(u))
end
function lambda(u) #微分方程式に出てくるラムダの項を計算する関数．ok
lamb =zeros(128,128)
for i in 1:39
for j in 1:128
lamb[i,j] = lam * (f[i,j]-u[i,j])
end #j
end #i
for i in 71:128
for j in 1:128
lamb[i,j] = lam * (f[i,j]-u[i,j])
end #j
end #i
lamb
end #function
#NLsolveを用いる
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
return nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini]).zero[1]
else
return nlsolve((vout,vin)->vout .= func(vin,params...), ini).zero
end
end
#nls(f,v,ini = u)
#非線形ソルバーで解くために方程式のスキームを書いた．
g(v)=(v -u)/dt + e * lap(lap(v)) - C1 * lap(v) + C2 * v -( lap((1/e) *4*(u^3) -6 * (u^2) + 2*u) + lambda(u) -C1 *(lap(u)) + C2*u)
nls(g,v, ini = u)
```
今書いてるやつ貼り付けますね
おｋ
これの最後で数値計算するんですがこけましたね
実行してみる。
どうでしたか？
うん、エラーがでな
実行結果どうなりましたか
```
g (generic function with 1 method)
```
だね。
定義しただけで終わってませんかそれ
たぶん
```
g(v)
```
を実行して終わった。
最後の行のnls()実行すれば計算がone step進められます．
ごめん
```
nls(g,v, ini = u)
```
忘れた。
あ，これだとvの定義忘れてます
dane
vの定義は
v = uです
おｋ
vはu_nに対してu_{n+1}の意味です
現状こんな感じ(めんどうだからgistに上げた)
ごめん、誤っていま消した。
もう一度発行する。
発行しなおした。
https://gist.github.com/muripoLife/500151545e6cb696a25314292d48a478
using NLsolve
が抜けていますね
https://gist.github.com/muripoLife/500151545e6cb696a25314292d48a478
いま、上げなおした。
エラーこれで同じ？
これです．同じです
じゃ、考えるか
エラーコードの
[12] (::var"#kw##nls")(::NamedTuple{(:ini,),Tuple{Array{Float64,2}}}, ::typeof(nls), ::Function, ::Array{Float64,2}) at ./none:0
にこけた原因が書かれてそうだと思いました（細かいところは読めてない．）
```
#NLsolveを用いる
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
return nlsolve((vout, vin)->vout[1] = func(vin[1], params...), [ini]).zero[1]
else
return nlsolve((vout, vin)->vout .= func(vin,params...), ini).zero
end
end
#nls(f,v,ini = u)
```
Numberって変数どこから出てきてる?
ボスのサイトにあったやつを丸写ししました
おｋ
http://www.cas.cmc.osaka-u.ac.jp/~paoon/misc/julia/post/nlsolve-package/
このサイトです
おｋ
これ結構時間かかる？
とりあえず、動いたみたいだけど、計算かえってこない
非線形方程式のソルバーなので時間はかかると思います．
いま動かしているコード
https://gist.github.com/muripoLife/500151545e6cb696a25314292d48a478
イテレーションごとにログ出ないの動いているかどうかわからんな
これはどこを書き換えましたか？
最期を
```
nls(g, ini=u)
```
にしただけ
更新されているかどうかわからない……
デバックしやすいように書きなおそ
確かに実行はされていそうですね
実行はされて、メモリーが5000mbかかっているところしかみてない。
いわゆるデバッグってどうするんですか
この場合、まずイテレーションごとに値プリントするとかかな？
今回は、gを高階関数で渡しているからかくの少し難しいかもね。
？？？
```
nls(g, ini=u)
```
gって関数でしょ。
関数を引数ように扱うのを高階関数というのよね。
合成関数とかがプログラミングで書きやすい的な。
なるほど
nlsの定義の中にうまく書き込んだらいいのかなと思うのですがどうでしょうか
それを試しているのだがうまく動かないね。
とりあえず長時間走らせてみます
ok
ワイもいろいろ書き換えてみるわ。
あとコード書いてるのjupyter?
そうです
30分やっても終わらないですね
おわらないね。
ちょいコードを考えておく。
```
using Images
using FFTW
using NLsolve
#種々のパラメーターのセッティング
dt = 1.0
e = 0.8
C1 = 300.0
C2 = 150000.0
lam = 50000.0
#テストするための画像を作る．
u = zeros(32,32) #真っ黒な画像
# make a view
img_a = @view u[1:32,8:13] #この行列を変化させると大元の行列も変えることができる．
img_a = fill!(img_a, 1.0)
img_b = @view u[1:32,20:24]
img_b = fill!(img_b, 1.0)
img_c = @view u[1:32,27:30]
img_c = fill!(img_c, 1.0)
img_d = @view u[16:20,1:32]
img_d = fill!(img_d, 0.5)
f = u #動かしたくない領域を参照するための行列
#uを更新していく
function nabla1(u)
fftu = fft(u)
fftu[1,:] = fftu[1,: ] * 0
for j in 1:15
fftu[j+1,: ] = (2 *pi *j *1im /32 ) * fftu[j+1 ,: ]
fftu[33-j ,: ] = (2 *pi *(-j) *1im /32 )*fftu[33-j ,: ]
end #j
real(ifft(fftu))
end
function nabla2(u)
fftu = fft(u)
fftu[: ,1 ] = fftu[: ,1 ] * 0
for k in 1:15
fftu[: ,k+1 ] = (2 *pi *k *1im /32 ) * fftu[: ,k+1 ]
fftu[: ,33-k ] = (2 *pi *(-k) *1im /32 )*fftu[: ,33-k ]
end #k
real(ifft(fftu))
end
function lap(u) #ラプラシアンの微分演算子．
nabla1(nabla1(u)) + nabla2(nabla2(u))
end
function lambda(u) #微分方程式に出てくるラムダの項を計算する関数．ok
lamb =zeros(32,32)
for i in 1:15
for j in 1:32
lamb[i,j] = lam * (f[i,j]-u[i,j])
end #j
end #i
for i in 21:32
for j in 1:32
lamb[i,j] = lam * (f[i,j]-u[i,j])
end #j
end #i
lamb
end #function
#NLsolveを用いる
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
return nlsolve((vout,vin)->vout[1]=func(vin[1],params...), [ini]).zero[1]
else
return nlsolve((vout,vin)->vout .= func(vin,params...), ini).zero
end
end
#nls(f,v,ini = u)
#非線形ソルバーで解くために方程式のスキームを書いた．vはuを1sステップ進めたもの
g(v)=(v -u)/dt + e * lap(lap(v)) - C1 * lap(v) + C2 * v -( lap((1/e) *4*(u^3) -6 * (u^2) + 2*u) + lambda(u) -C1 *(lap(u)) + C2*u)
for i in 1:10
@time u = nls(g, ini = u)
end
colorview(Gray,u)
```
ちゃんと動くようにはなりましたが，値にほぼ変化がないです．
Google Drive
Gist
Clipboard
Sign in with Wallet
