I have a bulky system of 26 differential equations with 26 variables. There is a trouble of solving such a system when I pass it to the function “ode!” as an array of symbolic equations. As a demonstration example of how it does not work, I write the following function
function ode!(du, initvals, system_of_eq, t)
eqs= vec(system_of_eq) #convert the system to a vector of equations
repl = vec(vcat([(z[j,i],initvals[j,i]) for j=1:2,i=1:3]...)) # i replace unknowns with the values of initvals
du = vec([eqs[k].subs(repl) for k=1:6]) # and substitute those values into the system
end
tspan = (0.0,30.0); #time span
init_val = [2.0 12.0 4.0; 1.0 1.0 1.0]; # initial conditions
z = [symbols("z_$j$i") for j = 1:2, i = 1:3]; #the array of unknowns
M = [z[1,2]+z[1,3]*z[1,2] z[1,1]+z[1,3]*z[1,2] z[2,3]-z[2,1]*z[2,2]; z[2,3]-z[2,1]*z[2,2] z[1,3] + z[1,1] z[2,3]*z[1,1] - z[1,3]] #the system of diff.equations dz^i/dt = f_i(z_1,...z_6)
Then I try to solve it:
prob = ODEProblem(od!,init_val,tspan,M)
sol = OrdinaryDiffEq.solve(prob,Tsit5(), reltol=1e-8, abstol=1e-8)
What I get is an array with the dimension (the size of initial values, number of steps) and not the solution.
Like this:
[:, :, 1] =
2.0 12.0 4.0
1.0 1.0 1.0
[:, :, 2] =
2.0 12.0 4.0
1.0 1.0 1.0
[:, :, 3] =
2.0 12.0 4.0
1.0 1.0 1.0
...
[:, :, 18] =
2.0 12.0 4.0
1.0 1.0 1.0
[:, :, 19] =
2.0 12.0 4.0
1.0 1.0 1.0
[:, :, 20] =
2.0 12.0 4.0
1.0 1.0 1.0
The same thing with the real system of 26 equations.
Are there any ideas why the function ‘solve’ does not work in this code? Another question is, how to change the step and the number of points which I need to find on the interval? By default, it searches for 20.
Any idea would be appreciated. Thank you.