 # ODE usage: passing parameters used for indexing

Hi all, I’m trying to solve a 2nd order (mass-spring-damper) ODE by splitting it into two first orders.

In the process of accumulating the force vector in my `accel()` function, I need to access 4 lists, a points list, a connectivity list, an edge list, and a resting-spring-length list. These aren’t actually parameters of the ODE function, but are necessary for accumulating the force vector. This is, maybe the problem, I’m not sure?

Here is the error I’m getting:

``````ERROR: LoadError: DimensionMismatch("expected axes([1.0, 0.0, 0.0], 1) = Base.OneTo(1)")
``````

Here is the code for the `accel()` function:

``````function accel(a, s, params, t)

PL, CL, EL, L, m, f, c = params

for el = 1:size(CL,2)

ea = (PL[:, CL[2,el]]+s[3*CL[2,el].-[2; 1; 0;]]) - (PL[:, CL[1,el]]+s[3*CL[1,el].-[2; 1; 0;]])

eb = (PL[:, CL[3,el]]+s[3*CL[3,el].-[2; 1; 0;]]) - (PL[:, CL[1,el]]+s[3*CL[1,el].-[2; 1; 0;]])

fe = (p/6)*cross(ea,eb)

ne = [3*CL[1,el].-[2; 1; 0] 3*CL[2,el].-[2;1;0] 3*CL[3,el].-[2; 1; 0]]

f[ne] = repeat(fe,3,1)

end

# Add spring & damping forces

for ele = 1:size(EL,2)

P = [ PL[:, EL[1,ele]]+s[3*EL[1,ele].-[2; 1; 0;]]   PL[:, EL[2,ele]]+s[3*EL[2,ele].-[2; 1; 0;]] ] # Coords of nodes belonging to this spring

V = [ s[(3*EL[1,ele].-[2; 1; 0;]).+96]     s[(3*EL[2,ele].-[2; 1; 0;]).+96] ] # Velocities of nodes belonging to this spring

dL = sqrt(sum((P[:,2] - P[:,1]).^2)) # Length of deformed spring (after t = 0)

disp = dL - L[ele] # Displacement from rest

relV = (V[:,1] - V[:,2]) # Relative velocity

n_hat = (P[:,2] - P[:,1])/dL # Directional unit vector

spr = L[ele]*disp # Spring force

damp = c[ele]*relV*n_hat # Damping force

ne = [3*EL[1,ele].-[2; 1; 0] 3*EL[2,ele].-[2; 1; 0]] # Indexing for each node within force vector

f[ne[:,1]] .+= -(spr + damp)*n_hat # For node1 of edge, f = k*(norm(pt2-pt1) - restLength)*(pt2-pt1)/norm(pt2-pt1)

f[ne[:,2]] .+= -(spr + damp)*-n_hat # Same for node2 of edge, except direction flipped, so that f = k*(norm(pt2-pt1) - restLength)*(pt1-pt2)/norm(pt2-pt1)

end

a[1:96] = s[1:96]

a[97:192] .= f ./ m

return a

end

``````

And my call:

``````params = (PL, CL, EL, L, m, f, c)

prob = ODEProblem(accel, s₀, (0.0,10.0), params)

sol = solve(prob)
``````

Is the error generated because I’m passing in those 4 lists, that don’t have the form of my output vector?

Any insight is appreciated, thank you!

Passing lists like that is fine, the parameters variable is just passed through the DifferentialEquations machinery unmodified so you can format it however you like.

That error looks like a broadcast between different size arrays, unfortunately I can’t guess where. Would you mind posting the full stack trace?

You could also try calling your `accel` method yourself with your initial conditions, that might produce an easier to understand stack trace.

``````du = similar(s₀)

accel(du, s₀, params, 0.0)
``````

Ok, yeah you’re right, it’s a broadcast error.

Thanks for the tip on calling `accel`, much easier to read the stack trace!