Adding parameters

We can also build up more complicated systems with multiple dependent variables and parameters as follows

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@parameters Dn, Dp
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs = [Dt(u(t, x)) ~ Dn * Dxx(u(t, x)) + u(t, x) * v(t, x),
    Dt(v(t, x)) ~ Dp * Dxx(v(t, x)) - u(t, x) * v(t, x)]
bcs = [u(0, x) ~ sin(pi * x / 2),
    v(0, x) ~ sin(pi * x / 2),
    u(t, 0) ~ 0.0, Dx(u(t, 1)) ~ 0.0,
    v(t, 0) ~ 0.0, Dx(v(t, 1)) ~ 0.0]

domains = [t ∈ Interval(0.0, 1.0),
    x ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem(
    eqs, bcs, domains, [t, x], [u(t, x), v(t, x)], [Dn => 0.5, Dp => 2])

discretization = MOLFiniteDifference([x => 0.1], t)

prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent

sol = solve(prob, Tsit5())

discrete_x = sol[x]
discrete_t = sol[t]

solu = sol[u(t, x)]
solv = sol[v(t, x)]

using Plots

anim = @animate for i in 1:length(t)
    p1 = plot(discrete_x, solu[i, :], label = "u, t=$(discrete_t[i])[1:9] ";
        legend = false, xlabel = "x", ylabel = "u", ylim = [0, 1])
    p2 = plot(discrete_x, solv[i, :], label = "v, t=$(discrete_t[i])";
        legend = false, xlabel = "x", ylabel = "v", ylim = [0, 1])
    plot(p1, p2)
end
gif(anim, "plot.gif", fps = 30)
Example block output

Remake with different parameter values

The system does not need to be re-discretized every time we want to plot with different parameters, the system can be remade with new parameters with remake. See the ModelingToolkit.jldocs for more ways to manipulate a prob post discretization.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@parameters Dn, Dp
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs = [Dt(u(t, x)) ~ Dn * Dxx(u(t, x)) + u(t, x) * v(t, x),
    Dt(v(t, x)) ~ Dp * Dxx(v(t, x)) - u(t, x) * v(t, x)]
bcs = [u(0, x) ~ sin(pi * x / 2),
    v(0, x) ~ sin(pi * x / 2),
    u(t, 0) ~ 0.0, Dx(u(t, 1)) ~ 0.0,
    v(t, 0) ~ 0.0, Dx(v(t, 1)) ~ 0.0]

domains = [t ∈ Interval(0.0, 1.0),
    x ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem(
    eqs, bcs, domains, [t, x], [u(t, x), v(t, x)], [Dn => 0.5, Dp => 2.0])

discretization = MOLFiniteDifference([x => 0.1], t)

prob = discretize(pdesys, discretization) # This gives an ODEProblem since it's time-dependent

sols = []
for (Dnval, Dpval) in zip(rand(10), rand(10))
    newprob = remake(prob, p = [Dnval, Dpval])
    push!(sols, solve(newprob, Tsit5()))
end

using Plots
for (j, sol) in enumerate(sols)
    discrete_x = sol[x]
    discrete_t = sol[t]
    solu = sol[u(t, x)]
    solv = sol[v(t, x)]
    anim = @animate for i in 1:length(discrete_t)
        p1 = plot(discrete_x, solu[i, :], label = "u, t=$(discrete_t[i])";
            legend = false, xlabel = "x", ylabel = "u", ylim = [0, 1])
        p2 = plot(discrete_x, solv[i, :], label = "v, t=$(discrete_t[i])";
            legend = false, xlabel = "x", ylabel = "v", ylim = [0, 1])
        plot(p1, p2)
    end
    gif(anim, "plot_$j.gif", fps = 10)
end
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
┌ Warning: : no method matching get_unit for arguments (Pair{Symbolics.Num, Float64},).
@ ModelingToolkit /cache/julia-buildkite-plugin/depots/01850c12-4ba3-484f-a49a-d388134cb7a1/packages/ModelingToolkit/Gpzyo/src/systems/validation.jl:154
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_1.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_2.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_3.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_4.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_5.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_6.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_7.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_8.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_9.gif
[ Info: Saved animation to /cache/build/tester-amdci4-9/julialang/methodoflines-dot-jl/docs/build/tutorials/plot_10.gif