I tried to implement the modulation manually with a for loop:
k = 1
I = [1.0, 1000.0, 10000.0, 100000.0, 1000000.0, 9000000.0, 9400000.0,
9450000.0, 9460000.0, 9470000.0, 9480000.0, 9490000.0,
9500000.0,
10000000.0, 80000000.0, 90000000.0, 100000000.0, 1000000000.0,
10000000000.0, 100000000000.0]
for i in I
println(k, " > ", i)
# inoculum
condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += i
cb1 = DiscreteCallback(condition1, affect1!)
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(growth!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])
# plot
clf()
myPlot = plot(soln.t, getindex.(soln.u, 1)+getindex.(soln.u, 2),
linestyle = "-", color = "blue")
myPlot = plot(soln.t, getindex.(soln.u, 3),
linestyle = "-", color = "green")
ax = gca()
ax.set_yscale("log")
ax.set_xlim([1, tmax])
ax.set_ylim([1e0, 1e10])
legend(["Bacteria (total)", "Phages"],
loc="lower right")
xlabel("Time (h)")
ylabel("Concentration (1/mL)")
Title = string("Inoculation of ", round(i; digits = 1), " units")
title(Title)
fig = string("/home/gigiux/Documents/model/plot/treat/treat_", k, ".png")
savefig(fig,
dpi = 300, format = "png", transparent = false)
k = k+1
end
and I got this breakpoint:
So by giving more than 9 490 000 phages, the bacteria go extinct. But I could proceed to check the between 9 490 000 and 9 485 000 and so forth. This approach is a very time-consuming way of proceeding.
Can I do the same (actually better) with SDE? Also because I would like to take into account also the time of inoculation…