Как исправить «LoadError: DimensionMismatch (« не может широковещательный массив иметь меньше измерений »)» - PullRequest
Новая
       77

Как исправить «LoadError: DimensionMismatch (« не может широковещательный массив иметь меньше измерений »)»

0 голосов
/ 11 февраля 2019

Я хотел бы решить следующие два связанных дифференциальных уравнения численно: 

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}

с определенными начальными условиями phi_0 (1-dim массив с N записями) и k_0 (2-dim массив с NxN записями)

Я попытался это сделать: Используя diffrentialEquations.js, построить матрицу начальных начальных условий u0 = hcat (Phi_0, k_0) (2-мерный массив, Nx (N + 1)) и каким-то образом определить, чтопервое уравнение применяется к первому столбцу (в моем коде [:, 1]), а второе уравнение применяется к другим столбцам (в моем коде [:, 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)

Запуск этих строк кода приводит к этой ошибке: 

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

Я новичок в Джулии, и я не уверен, что я направляюсь в правильном направлении с этим.Пожалуйста, помогите мне!

1 Ответ

0 голосов
/ 11 февраля 2019

Прежде всего, отредактируйте первый пакет, который Distributions, а не Distribution, мне потребовалось некоторое время, чтобы найти ошибку xD

Основная проблема - .= в вашем первомуравнение.Когда вы делаете это, вы не просто присваиваете новые значения массиву, вы делаете view.Я не могу точно объяснить вам, что такое представление, но я могу вам сказать, что при таком назначении левая и правая стороны должны иметь одинаковый тип.

Например: 

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 
...

Результатом является view, а не Вектор.С этим синтаксисом левая и правая стороны должны иметь один и тот же тип, чего не происходит в вашем примере.Обратите внимание, что типы rand(5) и rand(5,1) различны в Юлии: первый - Array{Float64,1}, а другой - Array{Float64,2}.В вашем коде d[:,1] - это Array{Float64,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) - это Array{Float64,2}, поэтому он не работает.У вас есть два варианта: изменить знак равенства на: 

du[:,1] = ...

Или: 

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]

Первый вариант - это просто базовое назначение, второй вариант использует способ view исоответствует типам обеих сторон.

...