Hello,
I’ve been having issues with my code related to the descretization step.
I tried taking the derivative of a fraction Dx(S(t,x)/(p-2L(t,x))) but I got this error, SO I decided to use the quotient rule to expand the derivative so that I was only taking the derivative of each function one at a time then multiplying them together accordingly.
I don’t really know how else to work around this error. Wondering if anyone has any comments?
I’ve checked my code 15 times for errors and I’m 90% sure it should work but I don’t know how to work around the derivative of the fraction part.
The function should read:
Here is my code:
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, DifferentialEquations, Plots
# Method of Manufactured Solutions: exact solution
# model parameters
le=0.25
Q = 0.002
p = 0.5
Dc = 0.0002
mu = 6.0
Yu = 0.587
ku = 0.2
X_inf = 1000.0
d = 0.5
a = 0.0
Y_lambda = 10000
eta = Y_lambda/X_inf
kd = 0.4
flux(S,L) = mu*(S)*L/((ku*(p-2*L)+S))
g(S,L) = mu*(S/(ku*(p-2*L)+S))
# 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*(((Dxx(S(t,x))*(p-2*L(t,x))+Dx(S(t,x))*(-2)*Dx(L(t,x)))-((-2)*Dxx(L(t,x))*S(t,x)+Dx(S(t,x))*(-2)*Dx(L(t,x))))*(p-2*L(t,x)) +4*Dx(L(t,x))*((p-2*L(t,x))*Dx(S(t,x)) - S(t,x)*(-2)*Dx(L(t,x))))/(p-2*L(t,x))^3 - 2*max(0,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*(((Dxx(W(t,x))*(p-2*L(t,x))+Dx(W(t,x))*(-2)*Dx(L(t,x)))-((-2)*Dxx(L(t,x))*W(t,x)+Dx(W(t,x))*(-2)*Dx(L(t,x))))*(p-2*L(t,x)) +4*Dx(L(t,x))*((p-2*L(t,x))*Dx(W(t,x)) - W(t,x)*(-2)*Dx(L(t,x))))/(p-2*L(t,x))^3 + 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*max(0,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,le)) ~ 0.0,
W(0, x) ~ 0.0,
W(t, 0) ~ 0.5,
Dx(W(t,le)) ~ 0.0,
L(0, x) ~ 0.0015]
# Space and time domains
domains = [t ∈ Interval(0.0, 10.0),
x ∈ Interval(0.0, le)]
# 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.0)
# 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 λ", label="t=$(discrete_t[i])")
end
plot(p1, p2, p3, legend=:false)
And the error message:
