Merge pull request #346 from mjsheikh/tutorial

Tutorial

Author: ChrisRackauckas

docs/make.jl

@@ -15,6 +15,9 @@ makedocs(
         "Symbolic Tutorials" => [
             "symbolic_tutorials/mol_heat.md"
         ],
+        "Operator Tutorials" => [
+            "operator_tutorials/kdv.md"
+        ],
         "Operators" => [
             "operators/operator_overview.md",
             "operators/derivative_operators.md",
@@ -26,8 +29,8 @@ makedocs(
         ],
         "Nonlinear Derivatives" => [
             "nonlinear_derivatives/nonlinear_diffusion.md"
-        ]
-    ]
+        ]  
+     ]
 )
 
 deploydocs(

docs/src/operator_tutorials/kdv.md

@@ -0,0 +1,50 @@
+# Solving KdV Solitons with Upwinding Operators
+
+The KdV equation is of the form `uₜ + αuuₓ + βuₓₓₓ = 0`. Here we'll use `α = 6`, `β = 1` for 
+simplicity of the true solution expression.
+
+### 1-Soliton solution using Upwind Difference
+
+The analytical expression for the single soliton case takes the form `u(x,t) = (c/2)/cosh²(√c * ξ/2)`.
+
+`c > 0` (wave speed) ; `ξ  = x - c*t` (moving coordinate)
+
+```julia
+using Test
+using DiffEqOperators, OrdinaryDiffEq, LinearAlgebra
+
+# Space domain and grids
+N = 21
+Δx = 1/(N-1)
+c = 1
+x = -10:Δx:10;
+
+# solution of the single forward moving wave
+ϕ(x,t) = (1/2)*sech.((x .- t)/2).^2 
+
+# Discretizing the PDE at t = 0
+u0 = ϕ(x,0);
+du = zeros(size(x)); 
+
+# Declaring the Upwind operator with winding = -1 since the wave travels from left to right 
+A = UpwindDifference{Float64}(1,3,Δx,length(x),-1);
+
+# Defining the ODE problem
+function KdV(du, u, p, t)
+	bc = GeneralBC([0,1,-6*ϕ(-10,t),0,-1],[0,1,-6*ϕ(10,t),0,-1],Δx,3) 
+	mul!(du,A,bc*u)
+end
+
+single_solition = ODEProblem(KdV, u0, (0.,5.));
+
+# Solving the ODE problem 
+soln = solve(single_solition,Tsit5(),abstol=1e-6,reltol=1e-6);
+
+# Plotting the results, comparing Analytical and Numerical solutions 
+using Plots
+plot(ϕ(x,0), title  = "Single forward moving wave", yaxis="u(x,t)", label = "t = 0.0 (Analytic)")
+plot!(ϕ(x,2), label = "t = 2.0 (Analytic)")
+plot!(soln(0.0), label = "t = 0.0 (Numerical)",ls = :dash)
+plot!(soln(2.0), label = "t = 2.0 (Numerical)",ls = :dash)
+```
+![solution_plot](https://user-images.githubusercontent.com/39168576/111033316-a702a180-8436-11eb-9b8e-749ad9206639.jpeg)

src/derivative_operators/derivative_operator.jl

@@ -26,34 +26,37 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv
     coeff_func              :: F
 end
 
-function nonlinear_diffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int,
+struct nonlinear_diffusion!{N} end
+struct nonlinear_diffusion{N} end
+
+function nonlinear_diffusion!{N}(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int,
     p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
     nknots::Int) where {T<:Real, N}
     #q is given by bc*u , u being the unknown function
     #p is given as some function of q , p being the diffusion coefficient
 
     @assert approx_order>1 "approximation_order must be greater than 1."
     if first_differential_order > 0 
-    du .= (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p)
+    du .= (CenteredDifference{N}(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference{N}(second_differential_order,approx_order,dx,nknots)*p)
     else 
-    du .= q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p)
+    du .= q[2:(nknots + 1)].*(CenteredDifference{N}(second_differential_order,approx_order,dx,nknots)*p)
     end
 
     for l = 1:(second_differential_order - 1)
-    du .= du .+ binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p)
+    du .= du .+ binomial(second_differential_order,l)*(CenteredDifference{N}(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference{N}(second_differential_order - l,approx_order,dx,nknots)*p)
     end
 
-    du .= du .+ (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]
+    du .= du .+ (CenteredDifference{N}(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]
 
 end
 
 # An out of place workaround for the mutating version
-function nonlinear_diffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int,
+function nonlinear_diffusion{N}(second_differential_order::Int, first_differential_order::Int, approx_order::Int,
     p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
     nknots::Int) where {T<:Real, N}
 
     du = similar(q,length(q) - 2)
-    return nonlinear_diffusion!(du,second_differential_order,first_differential_order,approx_order,p,q,dx,nknots)
+    return nonlinear_diffusion!{N}(du,second_differential_order,first_differential_order,approx_order,p,q,dx,nknots)
 end
 
 struct CenteredDifference{N} end
@@ -264,6 +267,8 @@ end
 
 CenteredDifference(args...) = CenteredDifference{1}(args...)
 UpwindDifference(args...) = UpwindDifference{1}(args...)
+nonlinear_diffusion(args...) = nonlinear_diffusion{1}(args...)
+nonlinear_diffusion!(args...) = nonlinear_diffusion!{1}(args...)
 use_winding(A::DerivativeOperator{T,N,Wind}) where {T,N,Wind} = Wind
 diff_axis(A::DerivativeOperator{T,N}) where {T,N} = N
 function ==(A1::DerivativeOperator, A2::DerivativeOperator)