Hello, I am currently trying to solve an ode where I have u(t) as a vector of length 1024, and I have du/dt = -u + AMf(u), where M is a matrix of 1024x1024. I specify some random u0, and solve the problem. I get a solution, but when I plot u(0) against u(100), as I am solving up to t = 100, I get the initial solution again. This happens even for simple problems, such as du = -u. Could someone help to explain what is going wrong? I have added the code below. I tried to make it self-contained, but if I made a mistake please let me know.
using DifferentialEquations
using Plots
pl = Plots
## Work to set up matrix
function W(x)
1/sqrt(Ď) .* exp.(-x.^2) .- 1/(sigma*sqrt(Ď)) .* exp.(-(x./sigma).^2)
end
n = 2^10; L = 10*Ď; h = (2*L)/n;
x = -L .+ (0:n-1)*h
M = zeros(n,n);
y = h .* W(x)
start = n/2 + 1;
for i = 1:n
M[i, :] .= circshift(y, start + i - 2)
end
#### Directly specifying initial conditions and parameters
Ď = 1.5
A = 1.7
Îź = 10.0
θ = 0.5
p = [Ď, A, Îź, θ]
u0 = 0.3*rand(n,1)
## Set up firing rate
function f(x)
1/(1 .+ exp.(-mu .* x + theta)) - 1/(1 .+ exp(theta))
end
# Setting up the problem
@. rhs(u,p,t) = -u + A*M*f(u)
# Setting up trial problem
# @. rhs(du, u, p, t) = -u
# Time steps
tspan = (0.0,100.0)
## Solving problem
prob = ODEProblem(rhs, u0, tspan, p)
sol = solve(prob, reltol = 1e-6, saveat = 0.001)
final = sol.u[end]
p1 = pl.plot(final)
p2 = pl.plot(u0)
pl.plot(p1,p2, layout = (1,2))
As @Rohan_Rajagopal said, your example code doesnât quite work. You can open a fresh REPL and execute the code yourself and fix all the errors about undefined variables and then update it here.
But also without running the code, Iâm a bit confused about the equations youâre trying to solve. M is a matrix, but then f(u) should be a vector and M * f(u) the regular matrix-vector product I guess?
Then your rhs(u, p, t) should probably look like (assuming A is just a number; note the missing dot between M and f(u)):
rhs(u, p, t) = - u .+ A .* M * f.(u)
(EDIT: I accidentally put a * instead of +)
and f(x) should have no dots, i.e. just map scalars to scalars.
The code should work now. I thought since I put @. in front of my rhs I do not need to add the dots for the multiplication. The equation is -u (n,1 vector) + a scalar * nxn matrix * nx1 vector. The product of the matrix and f gives a vector, so Im just adding these vectors together with the scalar multiplication as well. I will try what you wrote, thanks
Unfortunately the issue persists. I should also note that since u is a vector, shouldnât f have the dots, as its input is a vector? Or is this remedied by the f.(u)?
Yes, thatâs almost correct, but the issue is that @. adds the dot to every function call in the expression. But the matrix-vector product M * x has no dot, itâs just *.
Writing M .* x for a NxN matrix and a Nx1 vector will result in something different:
julia> M = rand(3,3)
3Ă3 Matrix{Float64}:
0.160868 0.997983 0.681036
0.834087 0.619388 0.646875
0.403048 0.658366 0.372333
julia> x = [1.0, 2.0, 3.0]
3-element Vector{Float64}:
1.0
2.0
3.0
julia> M * x
3-element Vector{Float64}:
4.199942527126877
4.013490057872009
2.8367784827159843
julia> M .* x
3Ă3 Matrix{Float64}:
0.160868 0.997983 0.681036
1.66817 1.23878 1.29375
1.20914 1.9751 1.117
This may or may not explain what you are seeing.
The latter. You can of course write the function f(u) such that it takes a vector and operates on it (with the dots as in your current example), but itâs common practice and in my opinion also easier to write and read the code if you write functions mapping scalars and then broadcast when you need it. The base functions also do it like that, e.g.
Okay thanks for help. I changed the code and checked that M*f.(u) was an Nx1 matrix, which it was. But unfortunately the issue still persists. I will check to see if I have made mistakes like this elsewhere. Out of curiosity, in the example you showed, it says M * x is a 3-element Vector. Mine says Matrix, and when I check the size it says Nx1. Is this a problem that it sees it as a matrix instead of a vector?
Yes, a Vector{Float64} and a 3Ă1 Matrix{Float64} are different objects in Julia, but usually both represent what we would just call âvectorâ in maths. They might behave similarly, but they also have some important differences in code.
E.g. * again computes the ânormalâ matrix-vector product and spits out another vector, whereas .* broadcasts the vector over all columns of the matrix:
julia> M = rand(3,3)
3Ă3 Matrix{Float64}:
0.81235 0.683495 0.260397
0.760464 0.447985 0.302931
0.354402 0.05534 0.79264
julia> x = [1.0, 2.0, 3.0]; y = reshape(x, (3, 1))
3Ă1 Matrix{Float64}:
1.0
2.0
3.0
julia> M * y
3Ă1 Matrix{Float64}:
2.960533034374361
2.5652254115093327
2.843002725837259
julia> M .* y
3Ă3 Matrix{Float64}:
0.81235 0.683495 0.260397
1.52093 0.895969 0.605861
1.06321 0.16602 2.37792
Note how the result of M .* y is again a matrix with each column multiplied (element-wise) with [1,2,3].
Regarding your original question,
If you ran the commented-out code from your post for du = -u then you will not see any change in your state because you are not updating du in your function. There are two different ways of specifying the ODE right-hand-side function and you mixed them up I believe. Here is an example of both versions side-by-side:
What the actual code does or should do is still not completely clear to me, as it doesnât run in a REPL (the parameters are sometimes called sigma and sometimes Ď and they are defined after they are used. I tried to fix the issues below (also removing all the dots I thought are unnecessary and using a plain Vector instead of a nx1 Matrix) and I donât get a constant solution for u(t). But I have no idea if Iâm solving the equations you intended to solve.
new code
using DifferentialEquations
using Plots
pl = Plots
## Work to set up matrix
Ď = 1.5
A = 1.7
Îź = 10.0
θ = 0.5
n = 2^10
L = 10*Ď
h = (2*L)/n;
p = [Ď, A, Îź, θ]
function W(x)
1/sqrt(Ď) * exp(-x^2) - 1/(Ď*sqrt(Ď)) * exp(-(x/Ď)^2)
end
x = -L .+ (0:n-1) .* h
M = zeros(n,n);
y = h .* W.(x)
start = n/2 + 1;
for i = 1:n
M[i, :] .= circshift(y, start + i - 2)
end
#### Directly specifying initial conditions and parameters
u0 = 0.3*rand(n)
## Set up firing rate
function f(x)
1/(1 + exp.(-Ο * x + θ)) - 1/(1 + exp(θ))
end
# Setting up the problem
rhs(u,p,t) = -u + A .* M * f.(u)
# Setting up trial problem
#### This function doesn't work as expected, since you would need to do `du .= -u`
# @. rhs(du, u, p, t) = -u
# Time steps
tspan = (0.0,1.0)
## Solving problem
prob = ODEProblem(rhs, u0, tspan, p)
sol = solve(prob, reltol = 1e-6, saveat = 0.1)
final = sol.u[end]
p1 = pl.plot(final)
p2 = pl.plot(u0)
pl.plot(p1,p2, layout = (1,2))
PS: The manual has a section about array arithmetic and broadcasting that might be helpful:
Thank you for the links. I will take a look. I am using JupyterNotebook, so maybe this is the issue with the REPL. I tried the code you wrote but unfortunately I still get u(end) = u(0). I also noticed that when I do M * f.(0) I still get a matrix instead of a vector. Not sure how to remedy this, as in the code you sent earlier M * x was a 3-element vector,
Dear all, thank you for the help. I logged in today and restarted my kernel in JupyterNotebook and now the code works. I have had this happen a couple times now, where I change nothing and code goes from either working to not or vice versa, but hopefully this does not happen again.
That is an inherent weakness of Jupyter notebooks⌠Have you tried Pluto.jl notebooks? Theyâll always keep your notebook in a consistent state. I cannot recommend them enough!
