Setting up differential equations

I’ve the following function

function fun(dx,x,p,t)
    dx[:,:] .= mat1*x + mat2*x
end

Here, I’d like to scale mat1*x + mat2*x by the values present in a vector scale.

In MATLAB I would do (1./scale).*(mat1*x + mat2*x). I’m not sure how element wise scaling is done in Julia for performing matrix vector multiplication.

mat1 is of size (10 x 10) , x : (10 x 1) scale : (10 x 1)

Suggestions?

@ptoche Could you please help me with this?

Example:


using DifferentialEquations, BenchmarkTools

mat1=[
       1    -2     1     0     0     0     0     0     0     0;
       0     1    -2     1     0     0     0     0     0     0;
       0     0     1    -2     1     0     0     0     0     0;
       0     0     0     1    -2     1     0     0     0     0;
       0     0     0     0     1    -2     1     0     0     0;
       0     0     0     0     0     1    -2     1     0     0;
       0     0     0     0     0     0     1    -2     1     0;
       0     0     0     0     0     0     0     1    -2     1;
       ];

mat2 = [
        1    -1     0     0     0     0     0     0     0     0;
        0     1    -1     0     0     0     0     0     0     0;
        0     0     1    -1     0     0     0     0     0     0;
        0     0     0     1    -1     0     0     0     0     0;
        0     0     0     0     1    -1     0     0     0     0;
        0     0     0     0     0     1    -1     0     0     0;
        0     0     0     0     0     0     1    -1     0     0;
        0     0     0     0     0     0     0     1    -1     0;
        ];

x0 = [1.0,0,0,0,0,0,0,0,0,0]
scale= [1 2 3 2 3 1 2 3];
saveat = 0:0.01:5

function fun(dx,x,p,t)
    dx[1,:] .= 0
    dx[2:9,:] .= (1/scale)*(mat1*x + mat2*x)
    dx[10,:] .= 2*(x[end-1] - x[end])
end

prob = ODEProblem(fun,x0,(0.0,5.0))
sys = modelingtoolkitize(prob)
fastprob = ODEProblem(sys,x0,(0.0,5.0),jac=true)

# Explicit RK Methods
@btime sol = solve(fastprob,Tsit5()) # 16.700 μs (245 allocations: 40.28 KiB)

doesn’t run successfully.

Should I do vecdot(1/scale, mat1*x+mat2*x) ?

Do you mean (mat1+mat2)*x./scale ? Try it in the REPL with 2 by 2 matrices and a simple x vector to see if that’s what you want, before throwing it into the ODEProblem. :thinking:

You can also play around with (one(10)./scale) and similar.

1 Like