Hi all,
I want to know why the code in CalculusWithJulia produce different solution for the same Initial Value problem:
y' = 0.2y(16-y)
y(0)=3
The plot is great anyway, since CairoMakie is a bit difficult to handle, this one really comes in handy.
so this is my problem, why this code:
# https://jverzani.github.io/CalculusWithJuliaNotes.jl/ODEs/odes.html
using SymPy, CalculusWithJulia, LaTeXStrings,Plots
@syms x
@syms y()
println("First-Order Linear Differential Equations with Julia")
println("y' = 0.2y(16-y) is the differential equation")
println("")
# we can use diff(u(x),x) but here, for visual simplicity, use the Differential operator,
D = Differential(x)
eqn = D(y)(x) ~ 0.2 * y(x) * (16 - y(x))
out = dsolve(eqn)
# The limits are confirmed by investigating the limits of the right-hand:
limit(rhs(out), x => oo), limit(rhs(out), x => -oo)
# For the initial value problem y(0)=3
eq = rhs(out) # just the right hand side
C1 = first(setdiff(free_symbols(eq), (x))) # fish out constant
c1 = solve(eq(x=>0) - 3, C1)
eq(C1 => c1[1])
sol = out(C1 => c1[1])
x0, y0 = 0, 3
ivp_sol = dsolve(eqn, y(x), ics=Dict(y(x0) => y0))
println("the particular solution for y(0) = 3: ")
sympy.pretty_print(ivp_sol)
produce different IVP solution than this code:
using SymPy, LaTeXStrings,Plots,PrettyTables
@syms x,y,c
@syms f()
println("First-Order Linear Differential Equations with Julia")
println("y' = 0.2y(16-y) is the differential equation")
println("")
p1(x) = -(3.2-0.2y)*x^0
g1(x) = 0
μ1(x) = exp(integrate(p1(x),x))
y1(x) = (1/μ1(x))*(integrate((μ1(x)*g1(x)),x) + c)
println("The general solution: ")
sympy.pretty_print(y1(x))
println("")
# Computing IVP with SymPy
diffeqf = Eq(f(x).diff(x) -(3.2 -0.2f(x))*x^0, g1(x)); string(diffeqf)
ivp_f = dsolve(diffeqf, f(x), ics = Dict(f(0) => 3))
println("the particular solution for velocity, f(0) = 3: ")
sympy.pretty_print(ivp_f)
I am able to plot it with this code:
using CalculusWithJulia, Plots, SymPy
@syms x y
x0, y0 = 0, 0
F(y, x) = 0.2y*(16-y)
plot(legend=false)
vectorfieldplot!((x,y) -> [1, F(y,x)], xlims=(x0, 10), ylims=(y0-1, y0+30))
g(x) = 16 - 13*exp(-0.2x)
plot!(g, linecolor=:blue, linewidth=2)
correct me if I am wrong.