Working on a matrix instead of vector

[gersonjferreira]

I'm working on a 2D diffusion equation with a single time derivative on the l.h.s. and x and y derivatives on the r.h.s:

d[ f(x,y,t) ]/dt = H f(x,y,t)

Without using ODE, I solve this with simply as

f(x,y,t+1) = f(x,y,t) + dt H f(x,y,t)

And it is more efficient to write f(x,y,t) as a matrix with x as lines as y as columns. So my function F(f) receives f as a matrix and returns another matrix out=H.f. In F(f) I work out the space derivatives using Fourier transforms.

This is not compatible with the ODE code, right? ODE requires f to be a vector.

Would it be difficult to generalize ODE to accept f(x,y,t) to be a matrix?

[acroy]

This should work. There is actually an example in test/Interface-tests.jl (
https://github.com/JuliaLang/ODE.jl/blob/master/test/interface-tests.jl).
Maybe you can have a look and try it. Just let us know if you have further
questions.

On Mi., 10. Feb. 2016 at 18:08, Gerson J. Ferreira [email protected]
wrote:

I'm working on a 2D diffusion equation with a single time derivative on
the l.h.s. and x and y derivatives on the r.h.s:

d[ f(x,y,t) ]/dt = H f(x,y,t)

Without using ODE, I solve this with simply as

f(x,y,t+1) = f(x,y,t) + dt H f(x,y,t)

And it is more efficient to write f(x,y,t) as a matrix with x as lines as
y as columns. So my function F(f) receives f as a matrix and returns
another matrix out=H.f. In F(f) I work out the space derivatives using
Fourier transforms.

This is not compatible with the ODE code, right? ODE requires f to be a
vector.

Would it be difficult to generalize ODE to accept f(x,y,t) to be a matrix?

—
Reply to this email directly or view it on GitHub
#86.

[mauro3]

You probably know this: often when working with ODE solvers, 2D (or ND) data is represented as a vector as far as the ODE solver is concerned. How f then treats the vector internally can be independent of that. For instance you can use reshape inside f to create a view in form of a 2D array.

[gersonjferreira]

Yes, reshape would work, but it makes a copy of the data, which could be slow. I was trying to avoid using reshape or rewriting the code in terms of a vector. It is just must simpler to evaluate the r.h.s. as it is now.

The example on the link (https://github.com/JuliaLang/ODE.jl/blob/master/test/interface-tests.jl) is great. Is the trick defining a custom type? In my code y is matrix, while in the example it is the custom type CompSol.

Here's a simple example that illustrate my question:

using ODE

h = [5.0 1.0; -1.0 6.0]/10.0;
dt = 1e-2;

y0 = [0.0; 1.0];
# y0 = [0.0 1.0; 1.0 1.0];

function rhs(t, y)
    return h*y
end

# solution via Euler
y=y0;
for t=0:dt:1
    y = y + dt*rhs(t,y)
end

# solution via RK45
tout, yout = ode45(rhs, y0, [0,1]);

# compare solutions
[y yout[end]]

If I go with the first definition of y0 (a vector) it all goes well, but if I choose the second definition of y0 (a matrix), it gives me

WARNING: [a] concatenation is deprecated; use collect(a) instead
in depwarn at deprecated.jl:73
in oldstyle_vcat_warning at ./abstractarray.jl:29
[inlined code] from abstractarray.jl:32
in oderk_adapt at /home/gerson/.julia/v0.5/ODE/src/runge_kutta.jl:220
[inlined code] from abstractarray.jl:33 (repeats 44 times)
in oderk_adapt at /home/gerson/.julia/v0.5/ODE/src/runge_kutta.jl:220
while loading /home/gerson/tmp/test.jl, in expression starting on line 18

[mauro3]

For normal arrays reshape does not copy:

julia> a = ones(4); b = reshape(a, 2,2); b[1,1]= 2; a[1]
2.0

but its copying semantics are unspecified. But otherwise you can also look at views, etc. But the scalar approach could be easiest. Actually, what happens if you just pass in a 2D array? Maybe that works as single dimensional indexing is defined.

[gersonjferreira]

Beautiful!! I didn't know reshape works like that! This will probably solve my problem. I cannot test right now, but I'll try soon! Thanks.

[mauro3]

Cool. Close this if it is resolved.

[gersonjferreira]

I'm not sure yet. My first attempt didn't work. But let me try again and I'll reply in a few days.

[mauro3]

Any luck?

[gersonjferreira]

Not really. The workaround was getting clumsy, so I decided to solve the differential equation with a different method (using the split-step method instead of RK).

But the problem seems to be on how the ODE package stores the output. In the code below I store the data at each time instant concatenating as an array of arrays. I believe this is similar to what ODE does, but somehow doesn't naturally works if the initial condition is a matrix.

using ODE

h = [5.0 1.0; -1.0 6.0]/10.0;
tspan=linspace(0,1,30);
dt=tspan[2]-tspan[1];

y0 = [0.0 1.0; 1.0 1.0];

function rhs(t, y)
    return h*y
end

y=y0;
yt=Any[y0];
for t=tspan[1:(end-1)]
    y = y + dt*rhs(t,y)
    yt = [yt Any[y]];
end

[mauro3]

Yes, you were right, it didn't work, except with ode4ms and ode5ms. So you can either us above PR until it gets merged or use ode*ms. Thanks!