Created
March 19, 2018 19:51
-
-
Save genkuroki/659798ab6f7a78d0d65ef21f34d116d8 to your computer and use it in GitHub Desktop.
FFTを用いた熱方程式やKdV方程式などを数値解法の in-place 版 jshtml no output
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "metadata": { | |
| "toc": true | |
| }, | |
| "cell_type": "markdown", | |
| "source": "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#FFTを用いた熱方程式やKdV方程式などを数値解法の-in-place-版\" data-toc-modified-id=\"FFTを用いた熱方程式やKdV方程式などを数値解法の-in-place-版-1\"><span class=\"toc-item-num\">1 </span>FFTを用いた熱方程式やKdV方程式などを数値解法の in-place 版</a></span><ul class=\"toc-item\"><li><span><a href=\"#諸定義\" data-toc-modified-id=\"諸定義-1.1\"><span class=\"toc-item-num\">1.1 </span>諸定義</a></span><ul class=\"toc-item\"><li><span><a href=\"#パッケージの読み込み\" data-toc-modified-id=\"パッケージの読み込み-1.1.1\"><span class=\"toc-item-num\">1.1.1 </span>パッケージの読み込み</a></span></li><li><span><a href=\"#数式を使った解説\" data-toc-modified-id=\"数式を使った解説-1.1.2\"><span class=\"toc-item-num\">1.1.2 </span>数式を使った解説</a></span></li><li><span><a href=\"#FFT_Data-タイプの定義\" data-toc-modified-id=\"FFT_Data-タイプの定義-1.1.3\"><span class=\"toc-item-num\">1.1.3 </span>FFT_Data タイプの定義</a></span></li><li><span><a href=\"#プロット用の函数\" data-toc-modified-id=\"プロット用の函数-1.1.4\"><span class=\"toc-item-num\">1.1.4 </span>プロット用の函数</a></span></li></ul></li><li><span><a href=\"#熱方程式で以上の定義のテスト\" data-toc-modified-id=\"熱方程式で以上の定義のテスト-1.2\"><span class=\"toc-item-num\">1.2 </span>熱方程式で以上の定義のテスト</a></span></li><li><span><a href=\"#KdV方程式の数値解法の最適化\" data-toc-modified-id=\"KdV方程式の数値解法の最適化-1.3\"><span class=\"toc-item-num\">1.3 </span>KdV方程式の数値解法の最適化</a></span><ul class=\"toc-item\"><li><span><a href=\"#プロット用函数\" data-toc-modified-id=\"プロット用函数-1.3.1\"><span class=\"toc-item-num\">1.3.1 </span>プロット用函数</a></span></li><li><span><a href=\"#FFTを-FFT-*-v,-FFT-\\-v-の形式で利用した場合\" data-toc-modified-id=\"FFTを-FFT-*-v,-FFT-\\-v-の形式で利用した場合-1.3.2\"><span class=\"toc-item-num\">1.3.2 </span>FFTを <code>FFT * v</code>, <code>FFT \\ v</code> の形式で利用した場合</a></span></li><li><span><a href=\"#FFTを-in-place-で利用した場合\" data-toc-modified-id=\"FFTを-in-place-で利用した場合-1.3.3\"><span class=\"toc-item-num\">1.3.3 </span>FFTを in-place で利用した場合</a></span></li><li><span><a href=\"#KdV方程式:-初期条件-sin\" data-toc-modified-id=\"KdV方程式:-初期条件-sin-1.3.4\"><span class=\"toc-item-num\">1.3.4 </span>KdV方程式: 初期条件 sin</a></span></li><li><span><a href=\"#ベンチマークテスト\" data-toc-modified-id=\"ベンチマークテスト-1.3.5\"><span class=\"toc-item-num\">1.3.5 </span>ベンチマークテスト</a></span></li><li><span><a href=\"#KdV方程式:-2-soliton解\" data-toc-modified-id=\"KdV方程式:-2-soliton解-1.3.6\"><span class=\"toc-item-num\">1.3.6 </span>KdV方程式: 2-soliton解</a></span></li></ul></li><li><span><a href=\"#1次元のSchroedinger方程式\" data-toc-modified-id=\"1次元のSchroedinger方程式-1.4\"><span class=\"toc-item-num\">1.4 </span>1次元のSchroedinger方程式</a></span><ul class=\"toc-item\"><li><span><a href=\"#Gaussian-packet\" data-toc-modified-id=\"Gaussian-packet-1.4.1\"><span class=\"toc-item-num\">1.4.1 </span>Gaussian packet</a></span></li><li><span><a href=\"#自由粒子\" data-toc-modified-id=\"自由粒子-1.4.2\"><span class=\"toc-item-num\">1.4.2 </span>自由粒子</a></span></li><li><span><a href=\"#低い壁\" data-toc-modified-id=\"低い壁-1.4.3\"><span class=\"toc-item-num\">1.4.3 </span>低い壁</a></span></li><li><span><a href=\"#高い壁\" data-toc-modified-id=\"高い壁-1.4.4\"><span class=\"toc-item-num\">1.4.4 </span>高い壁</a></span></li><li><span><a href=\"#調和振動子\" data-toc-modified-id=\"調和振動子-1.4.5\"><span class=\"toc-item-num\">1.4.5 </span>調和振動子</a></span></li><li><span><a href=\"#$V(x)-=--100-\\exp(-x^2/100)$\" data-toc-modified-id=\"$V(x)-=--100-\\exp(-x^2/100)$-1.4.6\"><span class=\"toc-item-num\">1.4.6 </span>$V(x) = -100 \\exp(-x^2/100)$</a></span></li></ul></li><li><span><a href=\"#Smith方程式\" data-toc-modified-id=\"Smith方程式-1.5\"><span class=\"toc-item-num\">1.5 </span>Smith方程式</a></span><ul class=\"toc-item\"><li><span><a href=\"#Smith方程式:-初期条件-sin\" data-toc-modified-id=\"Smith方程式:-初期条件-sin-1.5.1\"><span class=\"toc-item-num\">1.5.1 </span>Smith方程式: 初期条件 sin</a></span></li><li><span><a href=\"#Smith方程式:-初期条件-KdV-2-soliton\" data-toc-modified-id=\"Smith方程式:-初期条件-KdV-2-soliton-1.5.2\"><span class=\"toc-item-num\">1.5.2 </span>Smith方程式: 初期条件 KdV 2-soliton</a></span></li></ul></li></ul></li></ul></div>" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "# FFTを用いた熱方程式やKdV方程式などを数値解法の in-place 版\n\n黒木玄\n\n2018-01-23\n\n* Copyrigtht 2018 Gen Kuroki\n* License: MIT https://opensource.org/licenses/MIT\n\nこのノートブックは\n\n* http://nbviewer.jupyter.org/gist/genkuroki/14b05a2cfa172fc5ea351641d4bdda11\n* http://nbviewer.jupyter.org/gist/genkuroki/26928d4a1ae2e912a3850ed1b31e2941\n\nの続き. このノートブックでは in-place でFFTを使うことによってメモリ効率を最適化することを行いたい. 例としてKdV方程式を扱う." | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "@show Sys.CPU_CORES\nFFTW.set_num_threads(Sys.CPU_CORES)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## 諸定義" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### パッケージの読み込み" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "using PyPlot\n\ndefault_figsize = (6.4, 4.8) # このデフォルトサイズは大き過ぎるので\nrc(\"figure\", figsize=(3,2.4)) # このように小さくしておく.\n\nusing PyCall\nconst animation = pyimport(\"matplotlib.animation\")\n\nfunction displayfile(mimetype, file)\n open(file) do f\n base64text = base64encode(f)\n display(\"text/html\", \"\"\"<img src=\"data:$mimetype;base64,$base64text\">\"\"\")\n end\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### 数式を使った解説\n\n周期 $L=2\\pi K$ を持つ函数を\n\n$$\nf(x) = \n\\sum_{k=-N/2+1}^{N/2} (-1)^{k-1} a_k\\;e^{2\\pi i(k-1)x/L} =\n\\sum_{k=-N/2+1}^{N/2} (-1)^{k-1} a_k\\;e^{i(k-1)x/K}\n$$\n\nの形の有限Fourier級数でよく近似できていると仮定する. これの導函数は\n\n$$\nf'(x) = \\sum_{k=-N/2+1}^{N/2} (-1)^{k-1} \\frac{i(k-1)}{K}a_k\\;e^{i(k-1)x/K}\n$$\n\nになる. すなわち, $f(x)$ を微分する操作は, $-N/2< k \\leqq N/2$ のとき $a_k$ に $i(k-1)/K$ をかける操作に対応している. $-N/2< k \\leqq N/2$ のとき $a_k$ に $i(k-1)/K$ をかける操作は(丸め誤差)を除いて, 有限Fourier級数の正確な微分を与えることに注意せよ. この計算方法は小さな $h$ について $(f(x+h/2)-f(x-h/2))/h$ のような方法で微分を求める方法よりも誤差が小さい.\n\nさらに $x$ を\n\n$$\n\\Delta x = \\frac{L}{N} = \\frac{2\\pi K}{N}, \\quad\nx_j = (j-1)\\Delta x -\\frac{L}{2} = L\\left(\\frac{j-1}{N} - \\frac{1}{2}\\right)\n$$\n\nで離散化したとする. このとき,\n\n$$\nf(x_j) = \\sum_{k=-N/2-1}^{N/2} a_k\\;e^{i(k-1)(j-1)/N}.\n$$\n\nこれは $a_{k+N}=a_k$ という仮定のもとで\n\n$$\nf(x_j) = \\sum_{k=1}^N a_k\\; e^{i(k-1)(j-1)/N}\n$$\n\nと書き直される. これはJulia言語の離散Fourier変換のスタイルに一致している. この表示のもとで, $f(x)$ を微分する操作は, $1\\leqq k\\leqq N/2$ のとき $i(k-1)/N$ を $a_k$ にかけ, $N/2+1\\leqq k\\leqq N$ のとき $i(k-1-N)/N$ を $a_k$ にかける操作に対応している. `FFT_Data` における `D` を成分ごとにかける操作はまさにこの操作になっている." | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "?ifft", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### FFT_Data タイプの定義" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# right-open interval [a,b)\n#\nroi(a, b, N) = collect(linspace(a, b-(b-a)/N, N))\n\n# x axis (Assume the priodic boundary condition with period L)\n#\nx_axis(L, N) = roi(-L/2, L/2, N)\n\n# k axis (k = the wave number)\n#\nfunction k_axis(N)\n @assert iseven(N)\n Ndiv2 = div(N,2)\n vcat(roi(0,Ndiv2,Ndiv2), roi(-Ndiv2,0,Ndiv2))\nend\n\n# FFT Data\n#\nmutable struct FFT_Data{T<:Real, S}\n K::T\n L::T\n N::Int64\n x::Array{T,1}\n k::Array{T,1}\n D::Array{Complex{T},1}\n D2::Array{Complex{T},1}\n D3::Array{Complex{T},1}\n D4::Array{Complex{T},1}\n FFT::S\nend\n\n# FFT \\ (D.^m .* (FFT * u)) = (d/dx)^m u\n#\nfunction FFT_Data(K, N)\n L = 2π*K\n x = x_axis(L, N)\n k = k_axis(N)\n D = im .* k ./ K\n FFT = plan_fft(Array{Complex{Float64},1}(N))\n FFT_Data(Float64(K), L, N, x, k, D, D.^2, D.^3, D.^4, FFT)\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### プロット用の函数" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# nearly zero\n#\nfunction imagisnealyzero(z)\n maximum(abs, imag(z))/(maximum(abs, real(z)) + 1e-15) < 1e-3\nend\n\n# plot complex number array u on x\n#\nfunction plot_u(x, u)\n plot(x, real(u), label=\"Re\", lw=0.8)\n imagisnealyzero(u) || plot(x, imag(u), label=\"Im\", lw=0.8, ls=\"--\")\n grid(ls=\":\")\n xticks(fontsize=8)\n yticks(fontsize=8)\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## 熱方程式で以上の定義のテスト" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# Application - the solution of the heat equation u_t = u_{xx}\n#\n# u(t) = exp(t (∂/∂x)^2) u(0)\n#\nsol_of_heat_eq(o::FFT_Data, u0, t) = o.FFT \\ (exp.(t .* o.D2) .* (o.FFT * u0))\n\nN = 2^8\nK = 3\no = FFT_Data(K, N);\n\nf0(x) = exp(-4(x+4)^2) + exp(-2(x-3)^2)\nu0 = f0.(o.x)\nt = roi(0,2,10)\n@time u = [sol_of_heat_eq(o, u0, t) for t in t]\n\nfigure(figsize=(10, 3))\nfor j in 1:10\n subplot(2,5,j)\n plot_u(o.x, u[j])\n ylim(-0.05,1.05)\nend\ntight_layout()", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## KdV方程式の数値解法の最適化\n\n$$\nu_t = -u_{xxx}-3(u^2)_x\n$$" | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### プロット用函数" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "function plot_u(o, u, nrows, j, ymin, ymax)\n subplot(nrows, 5, j+1)\n plot_u(o.x, u)\n ylim(ymin, ymax)\n xticks(fontsize=8)\n yticks(fontsize=8)\n title(\"j = $j\", fontsize=8)\nend\n\nfunction plot_KdV(o::FFT_Data, u, t; thin=10, ymin=-7.0, ymax=15.0)\n nplots = div(length(t), thin)\n nrows = div(nplots+4, 5)\n\n figure(figsize=(10, 1.75nrows))\n \n j = 0\n plot_u(o, @view(u[:,j+1]), nrows, j, ymin, ymax)\n\n for j in 1:nplots-1\n plot_u(o, @view(u[:,thin*j+1]), nrows, j, ymin, ymax)\n end\n \n tight_layout()\nend\n\nfunction anim_KdV(figname, o::FFT_Data, u, t; thin=2, interval=100, ymin=-7.0, ymax=15.0)\n function plotframe(j)\n clf()\n plot(o.x, real(u[:, (j-1)*thin + 1]), lw=1.0)\n grid(ls=\":\")\n ylim(ymin, ymax)\n xticks(fontsize=8)\n yticks(fontsize=8)\n title(\"$(figname): j = $(j-1)\", fontsize=10)\n plot()\n end\n fig = figure(figsize=(4, 3))\n n = (length(t)-1)÷thin + 1\n frames = [1:n;]\n ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)\n anim_html = ani[:to_jshtml]()\n clf()\n display(\"text/html\", anim_html)\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### FFTを `FFT * v`, `FFT \\ v` の形式で利用した場合" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# FFT*v and FFT\\v version\n\nfunction update_KdV(o::FFT_Data, u, Δt, niters)\n v = o.FFT * u\n for iter in 1:niters\n v .= exp.(-Δt .* o.D3) .* v\n v .-= 3.*Δt .* o.D .*(o.FFT * (o.FFT \\ v).^2)\n end\n o.FFT \\ v\nend\n\nfunction solve_KdV(o::FFT_Data{T}, f0, tmax) where T\n Δt = 1/o.N^2\n skip = floor(Int, 0.005/Δt)\n t = 0:skip*Δt:tmax\n M = length(t)\n u = Array{Complex{T},2}(o.N, M)\n u[:,1] = Complex.(f0.(o.x))\n for j in 2:M\n u[:,j] = update_KdV(o, @view(u[:,j-1]), Δt, skip)\n end\n return real(u), t\nend\n\nN = 2^8\nK = 20/(2π)\no = FFT_Data(K, N);\n\nf0(x) = -5*sin(π*x/10)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.0\n@time u, t = solve_KdV(o, f0, tmax)\n\nsleep(0.1)\nplot_KdV(o, u, t)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### FFTを in-place で利用した場合" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "function update_KdV!(o::FFT_Data, u1, u0, Δt, niters)\n v = similar(u0)\n v_tmp = similar(u0)\n u_tmp = similar(u0)\n \n A_mul_B!(v, o.FFT, u0)\n for iter in 1:niters\n # v .= exp.(-Δt .* o.D3) .* v\n v .= exp.(-Δt .* o.D3) .* v\n \n # v .-= 3.*Δt .* o.D .*(o.FFT * (o.FFT \\ v).^2)\n A_ldiv_B!(u_tmp, o.FFT, v)\n u_tmp .*= u_tmp\n A_mul_B!(v_tmp, o.FFT, u_tmp)\n v .-= 3.*Δt .* o.D .* v_tmp\n end\n A_ldiv_B!(u1, o.FFT, v)\nend\n\nfunction solve_KdV_inplace(o::FFT_Data{T}, f0, tmax) where T\n Δt = 1/o.N^2\n skip = floor(Int, 0.005/Δt)\n t = 0:skip*Δt:tmax\n M = length(t)\n u = Array{Complex{T},2}(o.N, M)\n u[:,1] = Complex.(f0.(o.x))\n for j in 2:M\n update_KdV!(o, @view(u[:,j]), @view(u[:,j-1]), Δt, skip)\n end\n return real(u), t\nend\n\nN = 2^8\nK = 20/(2π)\no = FFT_Data(K, N);\n\nf0(x) = -5*sin(π*x/10)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.0\n@time u, t = solve_KdV_inplace(o, f0, tmax)\n\nsleep(0.1)\nplot_KdV(o, u, t)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### KdV方程式: 初期条件 sin" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"KdV sin\"\n@time anim_KdV(figname, o, u, t)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### ベンチマークテスト" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "using BenchmarkTools\n\nN = 2^8\nK = 20/(2π)\no = FFT_Data(K, N)\n\nf0(x) = -5*sin(π*x/10)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.4", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "@benchmark u, t = solve_KdV(o, f0, tmax)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "@benchmark u, t = solve_KdV_inplace(o, f0, tmax)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "以上のように in-place 版は非 in-place 版よりも計算時間を十数パーセント短く, in-place版の消費メモリは非 in-pace 版の50分の1である." | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### KdV方程式: 2-soliton解" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "N = 2^8\nK = 20/(2π)\no = FFT_Data(K, N);\n\nKdVsoliton(c, a, x) = c/2*(sech(sqrt(c)/2*(x-a)))^2\nf0(x) = KdVsoliton(16.0, -8.0, x) + KdVsoliton(4.0, -2.0, x)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.0\n@time u, t = solve_KdV_inplace(o, f0, tmax)\n\nsleep(0.1)\nplot_KdV(o, u, t, ymin=-0.5, ymax=10.0)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"KdV 2-soliton\"\n@time anim_KdV(figname, o, u, t, ymin=-0.5, ymax=10.0)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## 1次元のSchroedinger方程式" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# V はポテンシャル函数\n#\nmutable struct Schroedinger_Data{T, S, R}\n o::FFT_Data{T,S}\n V::R\n v::Array{T,1}\nend\n\n# FFT_Data o とポテンシャル函数 V から Schroedinger_Data を作成する函数\n#\nSchroedinger_Data(o, V) = Schroedinger_Data(o, V, V.(o.x))\n\nfunction update_Schroedinger!(s::Schroedinger_Data, u1, u0, Δt, skip)\n o = s.o\n U_tmp = similar(u0)\n u1 .= u0\n\n for iter in 1:skip\n # exp.(-im .* Δt .* s.v) .* (o.FFT\\(exp.(0.5im .* Δt .* o.D2).*(o.FFT * u)))\n A_mul_B!(U_tmp, o.FFT, u1)\n U_tmp .*= exp.(0.5im .* Δt .* o.D2)\n A_ldiv_B!(u1, o.FFT, U_tmp)\n u1 .*= exp.(-im .* Δt .* s.v)\n end\nend\n\nfunction solve_Schroedinger(s::Schroedinger_Data{T,S,R}, u0, tmax; Δt=2π/1000, skip=10) where {T,S,R}\n o = s.o\n t = 0:skip*Δt:tmax+10eps()\n M = length(t)\n u = Array{Complex{T},2}(o.N, M)\n u[:,1] = u0\n for j in 2:M\n update_Schroedinger!(s, @view(u[:,j]), @view(u[:,j-1]), Δt, skip)\n end\n return u, t\nend\n\nfunction plot_complexvalued(x, u, ymin, ymax)\n plot(o.x, real(u), lw=0.8)\n plot(o.x, imag(u), lw=0.8, ls=\"--\")\n ylim(ymin, ymax)\n grid(ls=\":\")\n xticks(fontsize=8)\n yticks(fontsize=8)\nend\n\nfunction plot_Schroedinger(s::Schroedinger_Data, u, t; thin=10, ymin=-1.1, ymax=1.1)\n nplots = (length(t)-1)÷thin + 1\n o = s.o\n nrows = div(nplots+3, 5)\n\n figure(figsize=(4, 1.75))\n subplot(121)\n plot_complexvalued(o.x, u[:,1], ymin, ymax)\n title(\"initial state\", fontsize=9)\n subplot(122)\n plot_u(o.x, s.v)\n title(\"potential\", fontsize=9)\n tight_layout()\n\n \n figure(figsize=(10, 1.75nrows))\n for j in 1:nplots-1\n subplot(nrows, 5, j)\n plot_complexvalued(o.x, u[:,(j-1)*thin+2], ymin, ymax)\n title(\"j = $j\", fontsize=9)\n end\n tight_layout()\nend\n\nfunction anim_Schroedinger(figname, s::Schroedinger_Data, u, t; thin=2, interval=100, ymin=-1.1, ymax=1.1, figtitle=figname)\n o = s.o\n function plotframe(j)\n clf()\n plot(o.x, real(u[:, (j-1)*thin+1]), lw=1.0)\n plot(o.x, imag(u[:, (j-1)*thin+1]), lw=1.0, ls=\"--\")\n grid(ls=\":\")\n ylim(ymin, ymax)\n xticks(fontsize=8)\n yticks(fontsize=8)\n title(\"$(figtitle): j = $(j-1)\", fontsize=10)\n plot()\n end\n fig = figure(figsize=(4, 3))\n n = (length(t)-1)÷thin + 1\n frames = [1:n;]\n ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)\n anim_html = ani[:to_jshtml]()\n clf()\n display(\"text/html\", anim_html)\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Gaussian packet" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "mutable struct GaussianPacket{T<:Real}\n x₀::T\n p₀::T\n σ::T\n Δx::T\nend\n\nfunction GaussianPacket(x₀, p₀, σ, o::FFT_Data)\n Δx = o.x[2] - o.x[1]\n GaussianPacket(x₀, p₀, σ, Δx)\nend\n\nfunction (f::GaussianPacket)(x)\n sqrt(f.Δx)/(π^(1/4)*sqrt(f.σ)) *\n exp(im*f.p₀*(x-f.x₀/2) - (x-f.x₀)^2/(2f.σ^2))\nend\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n\ng0 = GaussianPacket(-20.0, 5.0, 2.0, o)\nplot_u(o.x, g0.(o.x));", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### 自由粒子" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 自由粒子\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n@show typeof(o)\n\nV(x) = 0.0\ns = Schroedinger_Data(o, V)\n\ng0 = GaussianPacket(-20.0, 5.0, 2.0, o)\nu0 = g0.(o.x)\n\nT = 2.0\ntmax = 2π*T\n@time u, t = solve_Schroedinger(s, u0, tmax)\n\nsleep(0.1)\nplot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"free particle\"\n@time anim_Schroedinger(figname, s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### 低い壁" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 壁=(高さ10、幅10)、粒子=運動量5)\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n@show typeof(o)\n\nV(x) = ifelse(0 ≤ x ≤ 10, 10.0, 0.0)\ns = Schroedinger_Data(o, V)\n\ng0 = GaussianPacket(-20.0, 5.0, 1.0, o)\nu0 = g0.(o.x)\n\nT = 1.5\ntmax = 2π*T\n@time u, t = solve_Schroedinger(s, u0, tmax)\n\nsleep(0.1)\nplot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"low barrier\"\n@time anim_Schroedinger(figname, s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### 高い壁" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 壁=(高さ20、幅10)、粒子=運動量5)\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n@show typeof(o)\n\nV(x) = ifelse(0 ≤ x ≤ 10, 20.0, 0.0)\ns = Schroedinger_Data(o, V)\n\ng0 = GaussianPacket(-20.0, 5.0, 1.0, o)\nu0 = g0.(o.x)\n\nT = 1.5\ntmax = 2π*T\n@time u, t = solve_Schroedinger(s, u0, tmax)\n\nsleep(0.1)\nplot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"high barrier\"\n@time anim_Schroedinger(figname, s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### 調和振動子" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# 調和振動子 (周期 2π)\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n@show typeof(o)\n\nV(x) = x^2/2\ns = Schroedinger_Data(o, V)\n\ng0 = GaussianPacket(-20.0, 0.0, 0.5, o)\nu0 = g0.(o.x)\n\nT = 2.0\ntmax = 2π*T\n@time u, t = solve_Schroedinger(s, u0, tmax)\n\nsleep(0.1)\nplot_Schroedinger(s, u, t, ymin=-0.3, ymax=0.3)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"harmonic oscillator\"\n@time anim_Schroedinger(figname, s, u, t, ymin=-0.3, ymax=0.3, thin=1)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### $V(x) = -100 \\exp(-x^2/100)$" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "# V(x) = -20 exp(-x^2/25)\n\nN = 2^10\nK = 10\no = FFT_Data(K, N);\n@show typeof(o)\n\nV(x) = -100 * exp(-x^2/100)\ns = Schroedinger_Data(o, V)\n\ng0 = GaussianPacket(-15.0, 0.0, 1.0, o)\nu0 = g0.(o.x)\n\nT = 4.0\ntmax = 2π*T\n@time u, t = solve_Schroedinger(s, u0, tmax, skip=20)\n\nsleep(0.1)\nplot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"V(x)=-100exp(-x2／100)\"\nfigtitle = \"\\$V(x)=-100\\\\exp(-100x^2)\\$\"\n@time anim_Schroedinger(figname, s, u, t, ymin=-0.2, ymax=0.2, figtitle=figtitle)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "## Smith方程式\n\nKdV equation:\n$$\nu_t = -\\partial_x^3 u - 3\\partial_x(u^2).\n$$\n\nSmith's equation:\n$$\nu_t = 2a^{-2}\\left( \\sqrt{1-a^2\\partial_x^2} - 1 \\right)\\partial_x u - 3\\partial_x(u^2).\n$$\n\n以下では $a^2=0.15$ と仮定する." | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "function plot_Smith(o::FFT_Data, u, t; thin=10, ymin=-7.0, ymax=25.0)\n nplots = div(length(t), thin)\n nrows = div(nplots+4, 5)\n\n figure(figsize=(10, 1.75nrows))\n \n j = 0\n plot_u(o, @view(u[:,j+1]), nrows, j, ymin, ymax)\n\n for j in 1:nplots-1\n plot_u(o, @view(u[:,thin*j+1]), nrows, j, ymin, ymax)\n end\n \n tight_layout()\nend\n\nfunction anim_Smith(figname, o::FFT_Data, u, t; thin=2, interval=100, ymin=-7.0, ymax=25.0)\n function plotframe(j)\n clf()\n plot(o.x, real(u[:, (j-1)*thin + 1]), lw=1.0)\n grid(ls=\":\")\n ylim(ymin, ymax)\n xticks(fontsize=8)\n yticks(fontsize=8)\n title(\"$(figname): j = $(j-1)\", fontsize=10)\n plot()\n end\n fig = figure(figsize=(4, 3))\n n = (length(t)-1)÷thin + 1\n frames = [1:n;]\n ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)\n anim_html = ani[:to_jshtml]()\n clf()\n display(\"text/html\", anim_html)\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "function update_Smith!(o::FFT_Data, u1, u0, Δt, niters, a²)\n v = similar(u0)\n v_tmp = similar(u0)\n u_tmp = similar(u0)\n \n A_mul_B!(v, o.FFT, u0)\n for iter in 1:niters\n v .= exp.(2 ./ a² .* Δt .* (sqrt.(1 .- a² .* o.D2) .- 1).*o.D) .* v\n \n # v .-= 3.*Δt .* o.D .*(o.FFT * (o.FFT \\ v).^2)\n A_ldiv_B!(u_tmp, o.FFT, v)\n u_tmp .*= u_tmp\n A_mul_B!(v_tmp, o.FFT, u_tmp)\n v .-= 3.*Δt .* o.D .* v_tmp\n end\n A_ldiv_B!(u1, o.FFT, v)\nend\n\nfunction solve_Smith_inplace(o::FFT_Data{T}, f0, tmax; a²=0.15) where T\n Δt = 1/o.N^2\n skip = floor(Int, 0.005/Δt)\n t = 0:skip*Δt:tmax\n M = length(t)\n u = Array{Complex{T},2}(o.N, M)\n u[:,1] = Complex.(f0.(o.x))\n for j in 2:M\n update_Smith!(o, @view(u[:,j]), @view(u[:,j-1]), Δt, skip, a²)\n end\n return real(u), t\nend", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Smith方程式: 初期条件 sin" | |
| }, | |
| { | |
| "metadata": { | |
| "scrolled": false, | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "N = 2^9\nK = 20/(2π)\no = FFT_Data(K, N);\n\nf0(x) = -5*sin(π*x/10)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.0\n@time u, t = solve_Smith_inplace(o, f0, tmax)\n\nsleep(0.1)\nplot_Smith(o, u, t, thin=10)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"Smith sin\"\n@time anim_Smith(figname, o, u, t)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": {}, | |
| "cell_type": "markdown", | |
| "source": "### Smith方程式: 初期条件 KdV 2-soliton" | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "N = 2^9\nK = 20/(2π)\no = FFT_Data(K, N);\n\nKdVsoliton(c, a, x) = c/2*(sech(sqrt(c)/2*(x-a)))^2\nf0(x) = KdVsoliton(16.0, -8.0, x) + KdVsoliton(4.0, -2.0, x)\nΔt = 1/N^2\nskip = 10*floor(Int, 0.005/Δt)\n\ntmax = 1.0\n@time u, t = solve_Smith_inplace(o, f0, tmax)\n\nsleep(0.1)\nplot_Smith(o, u, t, ymin=-0.5, ymax=12.0)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "figname = \"Smith 2-soliton\"\n@time anim_Smith(figname, o, u, t, ymin=-0.5, ymax=12.0)", | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "metadata": { | |
| "trusted": true | |
| }, | |
| "cell_type": "code", | |
| "source": "", | |
| "execution_count": null, | |
| "outputs": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "name": "julia-0.6", | |
| "display_name": "Julia 0.6.2", | |
| "language": "julia" | |
| }, | |
| "toc": { | |
| "nav_menu": {}, | |
| "toc_window_display": false, | |
| "toc_section_display": true, | |
| "sideBar": true, | |
| "title_cell": "Table of Contents", | |
| "number_sections": true, | |
| "title_sidebar": "Contents", | |
| "skip_h1_title": false, | |
| "toc_position": { | |
| "width": "256px", | |
| "top": "150px", | |
| "left": "10px", | |
| "height": "calc(100% - 180px)" | |
| }, | |
| "toc_cell": true | |
| }, | |
| "_draft": { | |
| "nbviewer_url": "https://gist.github.com/d7328478c1188f876052fce859e91b40" | |
| }, | |
| "language_info": { | |
| "mimetype": "application/julia", | |
| "file_extension": ".jl", | |
| "version": "0.6.2", | |
| "name": "julia" | |
| }, | |
| "gist": { | |
| "id": "d7328478c1188f876052fce859e91b40", | |
| "data": { | |
| "description": "FFTを用いた熱方程式やKdV方程式などを数値解法の in-place 版 jshtml no output", | |
| "public": true | |
| } | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment