# Expansion of ODE equations

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 = dot(para01 * y, para02 * y) - y
dy = dot(para01 * y, exp.(para02 * y)) - y
dy = -dot(para01 * y, para02 * y) - y
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 += para01[ii] * y * para02[ii] * y
dy += para01[ii] * y * exp(para02[ii] * y)
dy += -para01[ii] * y * para02[ii] * y
end

dy -= y
dy -= y
dy -= y

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

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

2 Likes

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 += para01[ii] * y * para02[ii] * y
dy += para01[ii] * y * exp(para02[ii] * y)
dy += -para01[ii] * y * para02[ii] * y
end

dy -= y
dy -= y
dy -= y

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

Broadcast? `dy .= 0.0`

1 Like

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

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~

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

Can’t explain it better than the manual;)

1 Like

1 Like

Thanks for help

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

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

`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.

This post was temporarily hidden by the community for possibly being off-topic, inappropriate, or spammy.