DiffEqOperators and DifferentialEquations: solving a coupled system of PDEs

Specific Domains Numerics
4

You could try this, it works

using DiffEqOperators, OrdinaryDiffEq

# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
const bc = Dirichlet0BC(Float64)

function myModel!(du, u,p,t)
    du[1,:] .= Δ* bc * u[1,:]
    #du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
    du[2,:] .= Δ* bc * u[2,:]
end

# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (2, Nₓ))
prob = ODEProblem(myModel!, u0, tspan)
dt = 0.001
solve(prob, Tsit5(), dt = dt)

Since you have 2 PDEs, the initial condition should be u0 = rand(Float64, (2, Nₓ)) and you also need to match du & u in the function mymodel!

function myModel!(du, u,p,t)
    du[1,:] .= Δ* bc * u[1,:]
    #du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
    du[2,:] .= Δ* bc * u[2,:]
end
1 Like