Hello!
I’m working on solving PDEs in julia and producing plots of their solutions for different lengths of time.
I was having a lot of issues when using Julia v-1.9.0. When I increased the time in this specific setting past t=10, I had no issues with v-1.8.1 but the graphs fall apart when I use v-1.9.0 - see below
here is a copy of my code
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, DifferentialEquations, Plots, LinearAlgebra
# Method of Manufactured Solutions: exact solution
# model parameters
Q = 0.002
p = 0.5
Dc = 0.0002
mu = 6.0
Yu = 0.587
ku = 0.2
X_inf = 10000.0
d = 0.5
a = 0.5
Y_lambda = 0.587
eta = Y_lambda/X_inf
kd = 0.4
D = 0.0002
flux(S,L) = D*mu*((S/(p-2*L))/(ku+(S/(p-2*L))))*L
g(S,L) = mu*((S/(p-2*L))/(ku+(S/(p-2*L))))
# Parameters, variables, and derivatives
@parameters t x
@variables S(..), W(..), L(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
# 1D PDE and boundary conditions
eq = [Dt(S(t, x)) ~ (-Q*(p-2*L(t,x))*Dx(S(t,x)) + Q*S(t,x)*(-2)*Dx((L(t,x))))/(p-2*L(t,x))^2 + (Dc*(p-2*L(t,x))*Dxx(S(t,x)) - Dc*(-2)*Dxx(L(t,x))*(S(t,x)))/(p-2*L(t,x))^2 - 2*flux(S(t,x),L(t,x)) - (W(t,x)/Yu)*g(S(t,x),L(t,x)),
Dt(W(t, x)) ~ (-Q*(p-2*L(t,x))*Dx(W(t,x)) + Q*W(t,x)*(-2)*Dx((L(t,x))))/(p-2*L(t,x))^2 + (Dc*(p-2*L(t,x))*Dxx(W(t,x)) - Dc*(-2)*Dxx(L(t,x))*(W(t,x)))/(p-2*L(t,x))^2 + 2*X_inf*d*L(t,x) - 2*a*W(t,x) + g(S(t,x),L(t,x))*W(t,x),
Dt(L(t, x)) ~ eta*flux(S(t,x),L(t,x)) - (kd+d)L(t,x) + (a*W(t,x))/(X_inf)]
bcs = [S(0, x) ~ 5,
S(t, 0) ~ 5,
Dx(S(t,0)) ~ 0.0,
W(0, x) ~ 0.0,
W(t, 0) ~ 0.005,
Dx(W(t,0)) ~ 0.0,
L(0, x) ~ 0.00]
# Space and time domains
domains = [t ∈ Interval(0.0, 20.0),
x ∈ Interval(0.0, 1)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [S(t, x), W(t,x), L(t,x)])
# Method of lines discretization
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx], t)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
# Solve ODE problem
sol = solve(prob, Rosenbrock23(), saveat=2)
# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]
solL = sol[L(t,x)]
solC = sol[S(t,x)]./(p.-2 .*solL)
solU = sol[W(t,x)]./(p.-2 .*solL)
p1 = plot()
for i in 1:length(discrete_t)
plot!(discrete_x, solC[i, :], title="Numerical C", label="t=$(discrete_t[i])")
end
p2 = plot()
for i in 1:length(discrete_t)
plot!(discrete_x, solU[i,:], title="Numerical U", label="t=$(discrete_t[i])")
end
p3 = plot()
for i in 1:length(discrete_t)
plot!(discrete_x, solL[i, :], title="Numerical L", label="t=$(discrete_t[i])")
end
p = plot(p1, p2, p3, legend=:false)
Does anyone know why this is happening? Or how to make them work properly on v-1.9.0?
Thanks!
