Hello,
In one paragraph
I am trying to solve differential equations with discontinuities by solving them up to the discontinuity, change my system and input the old solution as initial conditions to a new problem. This gives me several output solutions and I need some way to combine them, either as solutions or by plotting them in succession.
In many paragraphs
I am working with chemical reaction networks and in my case the input to these are often of the forms of Steps/Pulses/Ramps. E.g. we have a system, suddenly a stimulus is turned on (a step), what is the response? Is there a difference if the response is turned on slowly, as opposed to a sudden increase (this would be a ramp). For simplicity I am not considering dirac pulses (maybe I should) but just a ramp up and then a ramp down shortly after,
I simulate my systems using Julias DifferentialEquations package (both ODEs and SDEs) and have understood that discontinuities like these should be handled with extra care. Initially I identified two ways to deal with them either using Callbacks (which seems complicated considering what I want to do, but maybe I am wrong?) or as my friend suggested: Solve the system until the discontinuity, use the final state of that solution as input to a new problem which is solved under the new conditions.
While the second initially seemed messy I think I could very easily make a small package which handles all of that to me. To it I would input my model, designate one or more parameters as changing and the package could do everything for me.
Due to the nature of my modelling DSL I thing this would also be advantageous as opposed to using callbacks, I might be wrong but I at least I think this solution would be simple enough and give me the desired results. (At least if the approach of solving the system multiple times with output of one solution becoming input to another is appropriate)
I have mostly figured out how I should do this, except for the final point were I have got stuck. The output of my program would be a set of solutions (sol1, sol2, …, solN). When I have finished my simulations I usually want to plot and investigate the result, but how do I do it with my disjoint solution? My two ideas would be to
- Preferably: Combine my solutions into a single solution, sol = (sol1 U sol2 U … U solN). I could then treat it as a normal problem output.
- Slightly more complicated: Plot solutions sol1, sol2, … sol3 in a single plot.
The first alternative I have gotten nowhere with. The second one slightly longer. If I have two solutions over the same interval sol1([0,1]) and sol2([0,1]) then I can plot them in the same plot using:
plot(sol1)
plot!(sol2)
However if they are over different intervals as in sol1([0,1]) and sol2([1,2]) (which is also the case in my intended application) then the plot window just jumps to the second interval once I write plot!(sol2) and I can no longer see my first solution. Surly there should be some way to solve this?
I am happy for all kind of input, whenever related to getting all plots displayed in the same window, or if I am doing something stupid in my entire approach.
If you are closer interested to know what I am doing here is more details:
I have a DSL which lets me input a reaction network like:
S = 0.1
p = [S]
rn = @reaction_network_new begin
S, X + Y--> XY
end S
were S is some constant, in this case defined to 0.1 I can then input rn and p to an ODEProblem which solves it for me. I want to have something which lets me do that, then changes the signal S to 1.0, input the old solution as initial conditions to a new ODEProblem and solves that, next I want to combine the two solutions somehow. (In this specific case I cannot do ramps, but if I can solve the problem for Steps I will move on to ramps later.)
