Expansion of ODE equations

New to Julia
,
1

New to julia, I have a question about how to expand the equations when solve the ODEs, here is a simple example
if the initial program is

function Eq(dy, y, p, t)
    para01, para02 = p
    dy[1] = dot(para01 * y[2], para02 * y[1]) - y[2]
    dy[2] = dot(para01 * y[1], exp.(para02 * y[2])) - y[3]
    dy[3] = -dot(para01 * y[3], para02 * y[1]) - y[1]
end

using DifferentialEquations, LinearAlgebra, Plots

para01 = [1.0, 2.0]
para02 = [3.0, 4.0]

tstart = 0
tend = 1e-2
tspan = (tstart, tend)
u0 = [0.0, 0.0, 1.0]
p = (para01, para02)
prob = ODEProblem(Eq, u0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol.t, sol[1, :])

Then I expand it to

function Eq(dy, y, p, t)
    para01, para02 = p

    for ii = 1:length(para01)
        dy[1] += para01[ii] * y[2] * para02[ii] * y[1]
        dy[2] += para01[ii] * y[1] * exp(para02[ii] * y[2])
        dy[3] += -para01[ii] * y[3] * para02[ii] * y[1]
    end

    dy[1] -= y[2]
    dy[2] -= y[3]
    dy[3] -= y[1]

end

using DifferentialEquations, LinearAlgebra, Plots

para01 = [1.0, 2.0]
para02 = [3.0, 4.0]

tstart = 0
tend = 1e-2
tspan = (tstart, tend)
u0 = [0.0, 0.0, 1.0]
p = (para01, para02)
prob = ODEProblem(Eq, u0, tspan, p)
sol = solve(prob, Tsit5())
plot!(sol.t, sol[1, :])

From my understanding, they should be the same, but the final result is different, and here is just a very simple case, for my real program, if I do this expansion, it can decrease the calculation time by a factor of 3, the problem is the result is different with and without expansion, any suggestions are welcome

2

You need to zero dy first, I don’t think there is any guarantee that it comes in as zero.

2 Likes
3

Thanks for reply, I know what you mean, but if I zero dy within the ODEs like this

function Eq(dy, y, p, t)
    para01, para02 = p
    dy=zeros(3)
    for ii = 1:length(para01)
        dy[1] += para01[ii] * y[2] * para02[ii] * y[1]
        dy[2] += para01[ii] * y[1] * exp(para02[ii] * y[2])
        dy[3] += -para01[ii] * y[3] * para02[ii] * y[1]
    end

    dy[1] -= y[2]
    dy[2] -= y[3]
    dy[3] -= y[1]

end

Then everything is wrong, if I zero dy out of this function, then it seems nothing happened, could you show me an example, how to do that

4

Broadcast? dy .= 0.0

1 Like
5

Yes, because dy=zeros(3) creates a new array but the solver uses still the old array. Then nothing happens.

6

Yes, if I put dy.=0.0 within the ODEs, it works, what’s the mechanism, I am still new to julia, will this occupy a lot of memory allocation, thanks~

7

https://docs.julialang.org/en/v1/manual/functions/#man-vectorized

Can’t explain it better than the manual;)

1 Like
8

Thanks for reply

1 Like
9

Thanks for help

10

That’s just equivalent to looping and doing dy[i] = 0, no memory allocations. fill!(dy,false) is also a good way to do it.

1 Like
11

Thanks for reply, both ways are working, and my ODEs now has 41 allocations without save everystep, I think it’s already good enough.
Finally I have one more question, if you have time, here is the example

a=rand(5)
b=rand(5)
c=rand(5)
d=rand(3)
e=rand(3)
dy=dot(a,b.*exp(c.+2))+dot(d,e)

initially without dot(d,e), I expand y=dot(a,b.*exp(c.+2)) with

for  i in 1:length(a)
    dy+=a[i]*b[i]*exp(c[i]+2)
end

but now with dot(d,e) can I use once for loop to put them together and get the same result as dy=dot(a,b.*exp(c.+2))+dot(d,e) , or maybe I can use the for loop twice

12

b.*exp(c.+2) is building an intermediate array which is then used in dot. Expanding it out is a better strategy for this kind of problem.