Wave equation on the square 3

Wave equation on the square 3.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "VERSION",
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 1,
          "data": {
            "text/plain": "v\"1.7.0-DEV.1129\""
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Wave equaiton: $\\displaystyle\n\\frac{\\partial^2 u}{\\partial t^2} = \n\\frac{\\partial^2 u}{\\partial x^2} + \n\\frac{\\partial^2 u}{\\partial y^2}\n$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "using DifferentialEquations\nusing Plots\nusing ProgressMeter\nusing Fmt\n\n# the wave equation on the square with free boundary condition\n\n#f!(du, v, u, p, t) = @. du = v\n\nfunction g!(dv, v, u, p, t)\n    h⁻² = p[1]\n    @. @views dv[2:end-1,2:end-1] = h⁻²*(\n        u[1:end-2,2:end-1] + u[3:end,2:end-1] + u[2:end-1,1:end-2] + u[2:end-1,3:end] - 4u[2:end-1,2:end-1])\n    @. @views dv[1  ,2:end-1] = h⁻²*(u[2,    2:end-1] + u[1,  1:end-2] + u[1,  3:end] - 3u[1,  2:end-1])\n    @. @views dv[end,2:end-1] = h⁻²*(u[end-1,2:end-1] + u[end,1:end-2] + u[end,3:end] - 3u[end,2:end-1])\n    @. @views dv[2:end-1,1  ] = h⁻²*(u[2:end-1,2    ] + u[1:end-2,1  ] + u[3:end,1  ] - 3u[2:end-1,1  ])\n    @. @views dv[2:end-1,end] = h⁻²*(u[2:end-1,end-1] + u[1:end-2,end] + u[3:end,end] - 3u[2:end-1,end])\n    dv[1  ,1  ] = h⁻²*(u[1  ,2    ] + u[2    ,1  ] - 2u[1  ,1  ]) \n    dv[end,1  ] = h⁻²*(u[end,2    ] + u[end-1,1  ] - 2u[end,1  ]) \n    dv[1  ,end] = h⁻²*(u[1  ,end-1] + u[2    ,end] - 2u[1  ,end]) \n    dv[end,end] = h⁻²*(u[end,end-1] + u[end-1,end] - 2u[end,end]) \nend\n\nx = y = range(-1, 1; length=401)\nh = step(x)\np = (1/h^2,)\n\nF(a, b, c, x, y) = exp(-((x-a)^2+(y-b)^2)/(2c^2))\nf(x, y) = F(0.7, 0, 0.01, x, y) + F(-0.4, 0.6, 0.01, x, y) + F(-0.4, -0.6, 0.01, x, y)\nu0 = f.(x', y)\nv0 = zero(u0)\ntspan = (0.0, 1.4)\n\n#prob = DynamicalODEProblem(g!, f!, v0, u0, tspan, p)\nprob = SecondOrderODEProblem(g!, v0, u0, tspan, p)\nsol = nothing; GC.gc()\n@time sol = solve(prob)\nsol = nothing; GC.gc()\n@time sol = solve(prob)\n\numin = minimum(minimum.(u.x[2] for u in sol.u))\numax = maximum(maximum.(u.x[2] for u in sol.u))\numin, umax",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": " 17.643820 seconds (31.82 M allocations: 10.026 GiB, 14.22% gc time)\n  9.177466 seconds (14.47 k allocations: 8.363 GiB, 21.04% gc time)\n",
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": "WARNING: Wrapping `Vararg` directly in UnionAll is deprecated (wrap the tuple instead).\n",
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "execution_count": 2,
          "data": {
            "text/plain": "(-0.2998091727151296, 1.0)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "L = 141\nts = range(tspan...; length=L) |> r -> [fill(r[1], 20); r; fill(r[end], 20)]\nprog = Progress(length(ts), 0)\n@gif for t in ts\n    heatmap(x, y, sol(t).x[2]';\n        clim=(-0.03, 0.03), colorbar=false, c=:bwr)\n    title!(f\"t = {$t:4.2f}\"; titlefontsize=12, aspectratio=1)\n    plot!(size=(420, 440), xlim=extrema(x), ylim=extrema(y))\n    next!(prog)\nend",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": "\u001b[32mProgress: 100%|█████████████████████████████████████████| Time: 0:01:22\u001b[39m\n┌ Info: Saved animation to \n│   fn = C:\\Users\\genkuroki\\OneDrive\\Math\\Math0052\\tmp.gif\n└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\iYDwd\\src\\animation.jl:104\n",
          "name": "stderr"
        },
        {
          "output_type": "execute_result",
          "execution_count": 3,
          "data": {
            "text/plain": "Plots.AnimatedGif(\"C:\\\\Users\\\\genkuroki\\\\OneDrive\\\\Math\\\\Math0052\\\\tmp.gif\")",