0

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!

hanf_lolo
  • 13
  • 4
  • The broadcast (dot) in `sum.([u[i,j+1] * sin(u[i,1] - u[j,1] + a) for i in 1:N, j in 1:N], dims=2)` seems suspicious. It will try to call `sum(::Float64, dims=2)` which should error as well. – carstenbauer Feb 11 '19 at 12:34
  • changed it, but the error message stayed the same. – hanf_lolo Feb 11 '19 at 13:35
  • Try to two extra lines: `temp = 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)` and `println("size of arrays: ", size(du), " and ", size(temp))` before the assignment to `du` and see what it prints. – DNF Feb 11 '19 at 13:50

1 Answers1

1

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.

tpdsantos
  • 423
  • 5
  • 6
  • It seems you are right, but I don't think it's a problem with `view`. The problem is that `broadcast` itself doesn't like assignment of a 2D source to a 1D destination. If you want broadcasted assignment, it is necessary to make sure that the output of `sum` is squeezed to 1D. – DNF Feb 11 '19 at 14:08
  • In that case `dropdims` is more efficient than `[:,1]`, but the simplest solution is to just drop the dot in the assignment. – DNF Feb 11 '19 at 14:19
  • didn't know about `dropdims` function, thanks. And after posting I remembered that you could just `vec()` it, probably more efficient as well. drop the dot is the simplest yeah, just put the second example in case OP wanted to mantain aesthetics of code. I would remove the dots of both equations, though. – tpdsantos Feb 11 '19 at 14:24