I want to know why the code in CalculusWithJulia produce different solution for the same Initial Value problem:
y' = 0.2y(16-y)
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) sol = out(C1 => c1) 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.