Solving a differential equation in a for loop

New to Julia
1

Hi Guys

I have done this before on Matlab however for some reason I can’t get it to run on Julia. I am getting the below error, I think it has something to do with regard to the for loop. Can someone please assist?

ERROR: ArgumentError: invalid index: 5.8e7 of type Float64

using DifferentialEquations
using Plots
function statespace!(dx, x, p,t)
    l1,l2,m,k,I,c = p
    dx[1] = x[2]
    dx[2] = (-k*x[1]*l1^2-c*x[2])/(m*l2^2+I)
end


x0 = [0.0, 0.025]  # initial conditions: rad and rad/s
tspan = (0,1)  # time span
k = [58*1000000.00,60*1000000] # N/m
m = 18000.00 # kg
I = 8500*5.15^2.00

for i in k
    p = (5.15,1.5,m,k[i],I,0.04*(k[i]+m))  # parameters: l1,l2,m,k,I,c
    prob = ODEProblem(statespace!, x0, tspan, p)
    sol = solve(prob)
    F_buffer = sol[1,:].*k[i]/1000
    plotly()
    plot(sol[1,:], xlabel="Time", ylabel="Displacement (m)")
    plot(F_buffer,xlabel = "Increment",ylabel = "Force in kN (Buffer)")
end
2

this is a vector of floats.
here you are trying to index an array using a float which is not allowed in Julia.

I think youre trying to do the following

for i in eachindex(k)

instead

note the the for loop loops over the elements of the array, not its indices

2 Likes
3

@ Salmon has already pointed out the issue.

I did some small changes and pasted result here for your given code.


using DifferentialEquations
using Plots
function statespace!(dx, x, p,t)
    l1,l2,m,k,I,c = p
    dx[1] = x[2]
    dx[2] = (-k*x[1]*l1^2-c*x[2])/(m*l2^2+I)
end


x0 = [0.0, 0.025]  # initial conditions: rad and rad/s
tspan = (0,1)  # time span
k = [58*1000000.00,60*1000000] # N/m
m = 18000.00 # kg
I = 8500*5.15^2.00

plt_1, plt_2 = plot(), plot()


for i in 1:length(k)
    p = (5.15,1.5,m,k[i],I,0.04*(k[i]+m))  # parameters: l1,l2,m,k,I,c
    prob = ODEProblem(statespace!, x0, tspan, p)
    sol = solve(prob)
    F_buffer = sol[1,:].*k[i]/1000

    plot!(plt_1, sol[1,:], xlabel="Time", ylabel="Displacement (m)")
    plot!(plt_2, F_buffer,xlabel = "Increment",ylabel = "Force in kN (Buffer)")
end

plot(plt_1, plt_2,  size = (1600, 850), thickness_scaling = 1.7, dpi = 300)

Result:

image
image846×448 41.1 KB

1 Like
4

for i in length(k) is equivalent to for i = 2, so it is only running the last version. I think for i in eachindex(k) is what you need.

I don’t use Plots.jl but you are plotting both lines on the same axis. The magnitude of one is 60_000 times the other which is making one of them appear to be unchanging. I think you want to put the two plots on different axes.

1 Like
5

Thank you for pointing that out I have edited the code. :slight_smile:

1 Like
6

Thank you all for the assistance