{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Plots,LaTeXStrings\n",
"default(markersize=3,linewidth=1.5)\n",
"using LinearAlgebra,DifferentialEquations\n",
"include(\"FNC.jl\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.1.1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the advection equation on a domain with periodic end conditions. Our approach is the method of lines."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"x,Dx,Dxx = FNC.diffper(300,[-4,4])\n",
"f = (u,c,t) -> -c*(Dx*u);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following initial condition isn't mathematically periodic, but the deviation is less than machine precision. We specify RK4 as the solver. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"u_init = @. 1 + exp(-3*x^2)\n",
"IVP = ODEProblem(f,u_init,(0.,3.),2.)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An animation shows the solution nicely. The bump moves with speed 2 to the right, reentering on the left as it exits to the right because of the periodic conditions. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/driscoll/Dropbox/books/fnc-extras/julia/advection.gif\n",
"└ @ Plots /Users/driscoll/.julia/packages/Plots/oiirH/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/driscoll/Dropbox/books/fnc-extras/julia/advection.gif\")"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"an = @animate for t = LinRange(0,3,120) \n",
" plot(x,sol(t),\n",
" xaxis=(L\"x\"),yaxis=([1,2],L\"u(x,t)\"), \n",
" title=\"Advection equation, t=$(round(t,digits=2))\",leg=:none )\n",
"end\n",
"gif(an,\"advection.gif\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.1.2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve for traffic flow using periodic boundary conditions. The following are parameters and a function relevant to defining the problem. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"rho_c = 1080; rho_m = 380; q_m = 10000;\n",
"Q0prime(rho) = q_m*4*rho_c^2*(rho_c-rho_m)*rho_m*(rho_m-rho)/(rho*(rho_c-2*rho_m) + rho_c*rho_m)^3;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we create a discretization on $m=800$ points."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"x,Dx,Dxx = FNC.diffper(800,[0,4]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we define the ODE resulting from the method of lines."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"ode = (rho,ep,t) -> -Q0prime.(rho).*(Dx*rho) + ep*(Dxx*rho);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our first initial condition has moderate density with a small bump."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rho_init = @. 400 + 10*exp(-20*(x-3)^2)\n",
"IVP = ODEProblem(ode,rho_init,(0.,1.),0.02)\n",
"sol = solve(IVP,alg_hints=[:stiff]);\n",
"\n",
"plot(x,[sol(t) for t=0:.2:1],label=[\"t=$t\" for t=0:.2:1],\n",
" xaxis=(L\"x\"),yaxis=(\"car density\"),title=\"Traffic flow\") "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The bump slowly moves backward on the roadway, spreading out and gradually fading away due to the presence of diffusion."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/driscoll/Dropbox/books/fnc-extras/julia/traffic1.gif\n",
"└ @ Plots /Users/driscoll/.julia/packages/Plots/oiirH/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/driscoll/Dropbox/books/fnc-extras/julia/traffic1.gif\")"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"an = @animate for t = LinRange(0,0.9,91) \n",
" plot(x,sol(t),\n",
" xaxis=(L\"x\"),yaxis=([400,410],\"density\"), \n",
" title=\"Traffic flow, t=$(round(t,digits=2))\",leg=:none )\n",
"end\n",
"gif(an,\"traffic1.gif\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we use an initial condition with a larger bump. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"rho_init = @. 400 + 80*exp(-16*(x-3)^2)\n",
"IVP = ODEProblem(ode,rho_init,(0.,0.5),0.02)\n",
"sol = solve(IVP,alg_hints=[:stiff]);"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/driscoll/Dropbox/books/fnc-extras/julia/traffic2.gif\n",
"└ @ Plots /Users/driscoll/.julia/packages/Plots/oiirH/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/driscoll/Dropbox/books/fnc-extras/julia/traffic2.gif\")"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"an = @animate for t = LinRange(0,0.5,101) \n",
" plot(x,sol(t),\n",
" xaxis=(L\"x\"),yaxis=([400,480],\"density\"), \n",
" title=\"A traffic jam, t=$(round(t,digits=2))\",leg=:none )\n",
"end\n",
"gif(an,\"traffic2.gif\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case the density bump travels backward along the road. It also steepens on the side facing the incoming traffic and decreases much more slowly on the other side. A motorist would experience this as an abrupt increase in density, followed by a much more gradual decrease in density and resulting gradual increase in speed. (You also see some transient, high-frequency oscillations. These are caused by instabilities, as we discuss in simpler situations later in this chapter.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.2.1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set up a test problem with velocity $c=2$ and periodic end conditions. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"x,Dx = FNC.diffper(400,[0 1])\n",
"uinit = @. exp(-80*(x-0.5)^2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this problem we use `RK4`, an explicit method. "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"ode = (u,c,t) -> -c*(Dx*u);\n",
"IVP = ODEProblem(ode,uinit,(0.,2.),2.)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = 2*(0:80)/80\n",
"u = hcat([sol(t) for t in t]...)\n",
"contour(x,t,u',color=:viridis,\n",
" xaxis=(L\"x\"),yaxis=(L\"t\"),title=\"Linear advection\",leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see the hump traveling rightward at constant speed, traversing the domain once for each integer multiple of $t=1/2$. We note the average time step that was chosen:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0035460992907801418"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"avgtau1 = sum(diff(sol.t))/(length(sol.t)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We cut $h$ by a factor of two and solve again."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"x,Dx = FNC.diffper(800,[0 1])\n",
"uinit = @. exp(-80*(x-0.5)^2);\n",
"IVP = ODEProblem(ode,uinit,(0.,2.),2.)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The CFL condition suggests that the time step should be cut by a factor of two also."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ratio = avgtau1 / avgtau2 = 1.9982269503546097\n"
]
},
{
"data": {
"text/plain": [
"1.9982269503546097"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"avgtau2 = sum(diff(sol.t))/(length(sol.t)-1)\n",
"@show ratio = avgtau1 / avgtau2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.2.2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set up advection over $[0,1]$ with velocity $c=-1$. This puts the right-side boundary in the upwind direction."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"n = 100\n",
"x,Dx = FNC.diffmat2(n,[0 1])\n",
"uinit = @. exp(-80*(x-0.5)^2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we try imposing $u=0$ at the right boundary, by appending that value to the end of the vector before multiplying by the differentiation matrix."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"trunc = u -> u[1:n]; extend = v -> [v;0];\n",
"ode = (v,c,t) -> -c*trunc(Dx*extend(v));\n",
"IVP = ODEProblem(ode,uinit[1:n],(0.,1.),-1)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(x[1:n],sol.t,sol[:,:]',\n",
" xlabel=L\"x\",ylabel=L\"t\",title=\"Inflow boundary condition\",leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data from the initial condition propagates out of the left edge. Because only zero is coming in from the upwind direction, the solution remains zero thereafter."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we try $u=0$ imposed at the left boundary. "
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"trunc = u -> u[2:n+1]; extend = v -> [0;v];\n",
"sol = solve(IVP,RK4());\n",
"plot(x[1:n],sol.t,sol[:,:]',\n",
" xlabel=L\"x\",ylabel=L\"t\",title=\"Outflow boundary condition\",leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Everything seems OK until the data begins to interact with the inappropriate boundary condition. The resulting \"reflection\" is entirely wrong for advection from right to left. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.3.1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $c=1$ we get the eigenvalues:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x,Dx = FNC.diffper(40,[0,1])\n",
"lambda = eigvals(Dx);\n",
"\n",
"scatter(real(lambda),imag(lambda),xlim=[-40,40],ylim=[-40,40],\n",
" title=\"Eigenvalues for pure advection\",leg=:none) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's choose a time step of $\\tau=0.1$ and compare to the stability regions of the Euler and backward Euler time steppers."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"zc = @.exp(1im*2pi*(0:360)/360); # points on |z|=1\n",
"\n",
"z = zc .- 1; # shift left by 1\n",
"plot(Shape(real(z),imag(z)),color=RGB(.8,.8,1),layout=(1,2),subplot=1)\n",
"scatter!(real(0.1*lambda),imag(0.1*lambda),subplot=1,\n",
" xlim=[-5,5],ylim=[-5,5],aspect_ratio=1,title=\"Euler\",leg=:none)\n",
"\n",
"z = zc .+ 1; # shift right by 1\n",
"plot!(Shape([-6,6,6,-6],[-6,-6,6,6]),color=RGB(.8,.8,1),subplot=2)\n",
"plot!(Shape(real(z),imag(z)),color=:white,subplot=2)\n",
"scatter!(real(0.1*lambda),imag(0.1*lambda),subplot=2,\n",
" xlim=[-5,5],ylim=[-5,5],aspect_ratio=1,title=\"Backward Euler\",leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the Euler case it's clear that *no* real value of $\\tau>0$ is going to make all (or even any) of the $\\tau\\lambda_j$ fit within the stability region. Hence Euler will never produce bounded solutions to this discretization of the advection equation. The A-stable backward Euler time stepping tells the exact opposite story; it will be absolutely stable regardless of $\\tau$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.3.2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The eigenvalues of advection--diffusion are near-imaginary for $\\epsilon\\approx 0$ and more negative-real for increasing values of $\\epsilon$."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot([],[],label=\"\",leg=:left,aspect_ratio=1)\n",
"x,Dx,Dxx = FNC.diffper(40,[0,1]);\n",
"tau = 0.1\n",
"for ep = [0.001 0.01 0.05]\n",
" lambda = eigvals(-Dx + ep*Dxx)\n",
" scatter!(real(tau*lambda),imag(tau*lambda),label=\"\\\\epsilon=$ep\")\n",
"end\n",
"title!(\"Eigenvalues for advection-diffusion\") "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.3.3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Deleting the first row and column places all the eigenvalues of the discretization into the left half of the complex plane. "
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x,Dx,Dxx = FNC.diffcheb(40,[0,1]);\n",
"A = -Dx[2:end,2:end]; # leave out first row and column\n",
"\n",
"lambda = eigvals(A)\n",
"scatter(real(lambda),imag(lambda),\n",
" xlim=[-300,100],title=\"Eigenvalues of advection with zero inflow\",leg=:none,aspect_ratio=1) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the rightmost eigenvalues have real part at most"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-4.931967035822995"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"maximum( real(lambda) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consequently all solutions decay exponentially to zero as $t\\to\\infty$. This matches the intuition of a flow with nothing at the inlet: eventually everything flows out of the domain. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.4.1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the wave equation (in Maxwell form) for speed $c=2$, with homogeneous Dirichlet conditions on the first variable. "
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"m = 200\n",
"x,Dx = FNC.diffmat2(m,[-1,1]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The boundary values of $u$ are given to be zero, so they are not unknowns in the ODEs we solve. Instead they are added or removed as necessary. "
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"trunc = u -> u[2:m];\n",
"extend = v -> [0;v;0];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function computes the time derivative of the system at interior points."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"dwdt = function(w,c,t)\n",
" u = extend(w[1:m-1])\n",
" z = w[m:2*m]\n",
" dudt = Dx*z\n",
" dzdt = c^2*(Dx*u)\n",
" return [ trunc(dudt); dzdt ]\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our initial condition is a single hump for $u$."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"u_init = @.exp(-100*(x)^2)\n",
"z_init = -u_init\n",
"w_init = [ trunc(u_init); z_init ]; "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because the wave equation is hyperbolic, we can use a nonstiff explicit solver."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"IVP = ODEProblem(dwdt,w_init,(0.,2.),2)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We extract the original $u$ and $z$ variables from the results, adding in the zeros at the boundaries for $u$."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"n = length(sol.t)-1\n",
"U = [ zeros(1,n+1); sol[1:m-1,:]; zeros(1,n+1) ];\n",
"Z = sol[m:2*m,:];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We plot the results for the original $u$ variable. "
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/driscoll/Dropbox/books/fnc-extras/julia/wave.gif\n",
"└ @ Plots /Users/driscoll/.julia/packages/Plots/oiirH/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/driscoll/Dropbox/books/fnc-extras/julia/wave.gif\")"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"an = @animate for (i,t) = enumerate(sol.t) \n",
" plot(x,U[:,i],layout=(1,2),subplot=1,\n",
" xaxis=(L\"x\"),yaxis=([-1,1],L\"u(x,t)\"), \n",
" title=\"Wave equation, t=$(round(t,digits=3))\",leg=:none )\n",
" plot!(x,sol.t,U',subplot=2,xlabel=L\"x\",ylabel=L\"t\",title=\"Space-time view\",leg=:none)\n",
"end\n",
"gif(an,\"wave.gif\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The original hump breaks into two pieces of different amplitudes. Each travels with speed $c=2$, and they pass through one another without interference. When a hump encounters a boundary, it is perfectly reflected, but with inverted shape. At time $t=2$ the exact solution looks just like the initial condition."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 12.4.2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now use a wave speed that is discontinuous at $x=0$. "
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"m = 120\n",
"x,Dx = FNC.diffcheb(m,[-1,1]);\n",
"c = @. 1 + (sign(x)+1)/2\n",
"trunc = u -> u[2:m];\n",
"extend = v -> [0;v;0];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function computes the time derivative of the method-of-lines system."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"dwdt = function(w,c,t)\n",
" u = extend(w[1:m-1])\n",
" z = w[m:2*m]\n",
" dudt = Dx*z\n",
" dzdt = c.^2 .* (Dx*u)\n",
" return [ trunc(dudt); dzdt ]\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set the initial conditions and solve using |ode45|."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"u_init = @.exp(-100*(x+0.5)^2);\n",
"z_init = -u_init;\n",
"w_init = [ trunc(u_init); z_init ]; \n",
"IVP = ODEProblem(dwdt,w_init,(0.,5.),c)\n",
"sol = solve(IVP,RK4());"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At each time in the following animation, we evaluate the discrete solution for $u$ and then extend it to the boundaries using the boundary conditions. "
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/driscoll/Dropbox/books/fnc-extras/julia/wavereflect.gif\n",
"└ @ Plots /Users/driscoll/.julia/packages/Plots/oiirH/src/animation.jl:90\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/driscoll/Dropbox/books/fnc-extras/julia/wavereflect.gif\")"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"an = @animate for t = LinRange(0,5,200)\n",
" plot(x,extend(sol(t,idxs=1:m-1)),\n",
" xaxis=(L\"x\"),yaxis=([-1,1],L\"u(x,t)\"), \n",
" title=\"Wave equation, variable speed\",label=\"t=$(round(t,digits=2))\" )\n",
"end\n",
"gif(an,\"wavereflect.gif\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each pass through the interface at $x=0$ generates a reflected and transmitted wave. By conservation of energy, these are both smaller in amplitude than the incoming bump. "
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.1.0",
"language": "julia",
"name": "julia-1.1"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.1.0"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}