CODEWITHSUNDEEP
indexingjuliaodedifferentialequations.jl
hanf_lolo
hanf_lolo

Reputation: 23

How to fix "LoadError: DimensionMismatch ("cannot broadcast array to have fewer dimensions")"

I'd like to solve the following two coupled differential equations numerically:

d/dt Phi_i = 1 - 1/N * \sum_{j=1}^N( k_{ij} sin(Phi_i - Phi_j + a) 
d/dt k_{ij} = - epsilon * (sin(Phi_i - Phi_j + b) + k_{ij}

with defined starting conditions phi_0 (1-dim array with N entries) and k_0 (2-dim array with NxN entries)

I tried this: Using DifferentialEquations.js, build a matrix of initial starting conditions u0 = hcat(Phi_0, k_0) (2-dim array, Nx(N+1)), and somehow define that the first equation applies to to first column (in my code [:,1]) , and the second equation applies to the other columns (in my code [:,2:N+1]).

using Distributions
using DifferentialEquations

N = 100
phi0 = rand(N)*2*pi
k0 = rand(Uniform(-1,1), N,N)

function dynamics(du, u, p, t)
    a = 0.3*pi 
    b = -0.53*pi 
    epsi = 0.01 

    du[:,1] .= 1 .- 1/N .* sum.([u[i,j+1] * sin(u[i,1] - u[j,1] + a) for i in 1:N, j in 1:N], dims=2)
    du[:,2:N+1] .= .- epsi .* [sin(u[i,1] - u[j,1] + b) + u[i,j+1] for i in 1:N, j in 1:N]
end

u0 = hcat(phi0, k0)
tspan = (0.0, 200.0)
prob = ODEProblem(dynamics, u0, tspan)
sol = solve(prob)

Running this lines of code result in this error: 

LoadError: DimensionMismatch ("cannot broadcast array to have fewer dimensions")in expression starting at line 47 (which is sol = solve(prob))

I'm new to Julia, and I'm not sure if im heading in the right direction with this. Please help me!

Upvotes: 0

Views: 1895

Answers (1)

tpdsantos
tpdsantos

Reputation: 443

First of all, edit the first package, which is Distributions and not Distribution, it took me a while to find the error xD

The main problem is the .= in your first equation. When you do that, you don't just assign new values to an array, you're making a view. I cannot explain you exactly what is a view, but what I can tell you is that, when you this kind of assign, the left and right side must have the same type.

For example:

N = 100
u = rand(N,N+1)
du = rand(N,N+1)

julia> u[:,1] .= du[:,1]
100-element view(::Array{Float64,2}, :, 1) with eltype Float64:
 0.2948248997313967 
 0.2152933893895821 
 0.09114453738716022
 0.35018616658607926
 0.7788869975259098 
 0.2833659299216609 
 0.9093344091412392 
...

The result is a view and not a Vector. With this syntax, left and right sides must have same type, and that does not happen in your example. Note that the types of rand(5) and rand(5,1) are different in Julia: the first is an Array{Float64,1} and the other is Array{Float64,2}. In your code, d[:,1] is an Array{Float64,1} but 1 .- 1/N .* sum.([u[i,j+1] * sin(u[i,1] - u[j,1] + a) for i in 1:N, j in 1:N], dims=2) is an Array{Float64,2}, that's why it doesn't work. You have two choices, change the equal sign for:

du[:,1] = ...

Or:

du[:,1] .= 1 .- 1/N .* sum.([u[i,j+1] * sin(u[i,1] - u[j,1] + a) for i in 1:N, j in 1:N], dims=2)[:,1]

The first choice is just a basic assign, the second choice uses the view way and matches the types of both sides.

Upvotes: 1

Related Questions