Dimension mismatch error

General Usage Performance
1

Hello everyone,

I’m working on a project in Julia to predict reaction mechanisms for a given set of inputs using kinetic modeling. I’m fairly new to Julia, and everything seems to be running well except for a dimension mismatch error in the plotting section.

The Issue

The problem arises because I’m using the diff function to calculate the reaction rate (i.e., the change in mass over time). This results in an array with a length of ( n - 1 ), whereas the temperature array, derived from the solution time points, has a length of ( n ). I attempted to adjust by truncating the temperature array (Temp[1:end-1]) to match the length of the diff output, but the error persists.

Code Snippet

Here’s the relevant part of my code and issue:

# Adjust Temp to match the length of `diff(mass_r1)` and `diff(time)`
plot(Temp[1:end-1] .- 273.15, -diff(mass_r1) ./ diff(time), label="R1", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r2) ./ diff(time), label="R2", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r3) ./ diff(time), label="R3", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r4) ./ diff(time), label="R4", linewidth=1.5)

DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 7 and 5

Stacktrace:
 [1] _bcs1
   @ .\broadcast.jl:529 [inlined]
 [2] _bcs
   @ .\broadcast.jl:523 [inlined]
 [3] broadcast_shape
   @ .\broadcast.jl:517 [inlined]
 [4] combine_axes
   @ .\broadcast.jl:512 [inlined]
 [5] instantiate
   @ .\broadcast.jl:294 [inlined]
 [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(/), Tuple{Vector{Float64}, Vector{Float64}}})
   @ Base.Broadcast .\broadcast.jl:873
 [7] top-level scope
   @ In[7]:3

My Question

What would be the best way to handle the dimension mismatch here? Is there a more effective or Julia-friendly approach to ensure that the arrays match up correctly? I’d like to know if there’s a way to plot reaction rates calculated with diff against temperature in Julia.

Thanks in advance for any suggestions!

2

It’s hard to help you unless you post runnable code, with example data for your arrays. It may be that the source of the problem is not diff.

3

Hello. here is the whole code snippet I used.

using DifferentialEquations
using LinearAlgebra
using DelimitedFiles
using Plots

const R = 8.314  # Universal gas constant in J/(mol*K)

# Define the function to calculate temperature over time with a constant heating rate
function getsampletemp(t, T0, beta)
    T = clamp.((T0 + beta / 60 * t), 0, 873.15)   # K/min to K/s
    return T
end

# Define the CRNN function without w_out
function crnn!(du, u, p, t)
    T = getsampletemp(t, p.T0, p.beta)

    # Extract kinetic parameters from p
    a = p.w_in[1]        # stoichiometric coefficient/reaction order
    b = p.w_in[2]        # non-temperature dependence (T^b term)
    Ea = p.w_in[3]       # activation energy (in J/mol)
    A = exp(p.w_b[1])    # pre-exponential factor

    # Calculate reaction rate constant k using the modified Arrhenius equation
    k = A * T^b * exp(-Ea / (R * T))

    # Calculate rates for reactant and product
    du[1] = -a * k * u[1]^a    # change in reactant concentration
    du[2] = a * k * u[1]^a     # change in product concentration
end

# Define initial parameters
mass = 5.0                   # Initial mass in mg
T0 = 290.0                   # Initial temperature in K
beta = 2.5 / 60              # Heating rate in K/s

# Define kinetic parameters for each reaction
# Reaction 1
w_in1 = [0.717, 0.26, 133270.0]   # [stoichiometric coefficient, non-temp dependence, activation energy]
w_b1 = [20.83]                     # pre-exponential factor

# Reaction 2
w_in2 = [0.24, 0.0, 189780.0]
w_b2 = [13.13]

# Reaction 3
w_in3 = [1.0, 0.13, 157620.0]
w_b3 = [15.36]

# Reaction 4
w_in4 = [0.68, 0.33, 148600.0]
w_b4 = [22.43]

# Initial concentrations
u0 = [mass, 0.0]   # [initial reactant mass, initial product mass]

# Define problem and parameters for each reaction
param1 = (w_in=w_in1, w_b=w_b1, T0=T0, beta=beta)
param2 = (w_in=w_in2, w_b=w_b2, T0=T0, beta=beta)
param3 = (w_in=w_in3, w_b=w_b3, T0=T0, beta=beta)
param4 = (w_in=w_in4, w_b=w_b4, T0=T0, beta=beta)

tspan = (0.0, 1.0)  # Define time span for integration (seconds)

# Define and solve the ODE problem for each reaction
prob1 = ODEProblem(crnn!, u0, tspan, param1)
sol1 = solve(prob1, Tsit5())

prob2 = ODEProblem(crnn!, u0, tspan, param2)
sol2 = solve(prob2, Tsit5())

prob3 = ODEProblem(crnn!, u0, tspan, param3)
sol3 = solve(prob3, Tsit5())

prob4 = ODEProblem(crnn!, u0, tspan, param4)
sol4 = solve(prob4, Tsit5())


# Extract data for plotting
time = sol1.t
Temp = getsampletemp.(time, T0, beta)  # Calculate temperature for each time point

mass_r1 = sol1.u |> x -> [xi[1] for xi in x]  # Reactant mass for Reaction 1
mass_p1 = sol1.u |> x -> [xi[2] for xi in x]  # Product mass for Reaction 1

mass_r2 = sol2.u |> x -> [xi[1] for xi in x]
mass_p2 = sol2.u |> x -> [xi[2] for xi in x]

mass_r3 = sol3.u |> x -> [xi[1] for xi in x]
mass_p3 = sol3.u |> x -> [xi[2] for xi in x]

mass_r4 = sol4.u |> x -> [xi[1] for xi in x]
mass_p4 = sol4.u |> x -> [xi[2] for xi in x]

# Adjust Temp to match the length of `diff(mass_r1)` and `diff(time)`
plot(Temp[1:end-1] .- 273.15, -diff(mass_r1) ./ diff(time), label="R1", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r2) ./ diff(time), label="R2", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r3) ./ diff(time), label="R3", linewidth=1.5)
plot!(Temp[1:end-1] .- 273.15, -diff(mass_r4) ./ diff(time), label="R4", linewidth=1.5)

xlabel!("Temperature [°C]")
ylabel!("-dm/dt [mg/s]")
legend!(location=:northeast)

I hope that helps.

4

It’s not ideal, because it requires people to install the whole DifferentialEquations stack (≈ 200 packages) and run your whole simulation just to see why one line of code isn’t working.

It would be better to just post sample data Temp = [...]; mass_r1 = [...]; time = [...].

5

I didn’t run the code, but I’m guessing that mass_r2, mass_r3, and mass_r4 all have different lengths compared to mass_r1 (which has the same length as time and temp.

You can evaluate an ODESolution at given a given timepoint by calling e.g., sol1(t). Thus, you could construct

mass_r2 = [sol2(ti)[1] for ti in time]

and make the other objects in a similar manner. The odesolver will not in general take the same timesteps when you change the input in some way (in your case the parameter set)

2 Likes
6

Right, mass_r2 is a different length:

julia> length.((time, Temp, mass_r1))
(6, 6, 6)

julia> length.((mass_r2, mass_r3, mass_r4))
(8, 6, 6)
7

On a separate note: using a finite difference to calculate the reaction rate is not ideal — you are throwing away a lot of the accuracy of the solution.

You can get the “exact” time derivative (or at least, as accurate as your solution) just by plugging the solution into your right-hand-side function crnn!. But even more simply you can do e.g.

dmass_r1_dt = getindex.(sol1.(time, Val{1}), 1) # Val{1} = first derivative

or replace the getindex call with

dmass_r1_dt = sol1.(time, Val{1}, idxs=1)

or

dmass_r1_dt = sol1(time, Val{1}, idxs=1).u

See also DifferentialEquations and derivatives of solutions? and “Interpolations and Calculating Derivatives” in the manual.

2 Likes
8

Thank you for the help. The provided code worked like wonders.

9

It was my first time asking in a community. I did not take that in consideration. I shall be careful from now on. Thank you

10

Thank you so much for the insight. I shall work on it. As I mentioned already, I am relatively new to both coding and that too in julia. And that applies same for using this community. If I have an additional question on a different section related to the same code, should I continue here or make a new thread?

11

Generally, we have one thread per problem here :slight_smile:
So if it’s (likely) a different problem, then make a new thread. If it requires a new explanation, then it certainly is a different problem.

Use this thread, if you want to ask more about the suggested solution or they don’t solve your problem.