Differential Equation problem with Index

Specific Domains Modelling & Simulations
1

Hello, I have a model which I basically need to solve through the calculation of many differential equations by examining the surrounding values and such.

I believe the equation is correctly being solved. However, there is one issue it has been having which effects future values.
Basically, for an the result of an equation like this, the value of 2.066x10^23 should be modified to have the sum of every non-negative value in the triangle shape subtracted from it.

 -1.0       -1.0           -1.0            -1.0          0.0
 -1.0       -1.0           -1.0             9.89029e-11  1.90176e-17
 -1.0       -1.0         1681.92            0.00290038   1.72268e-10
  6.022e23   2.86406e17     4.22592e10  12482.0          0.000225923

IE: 2.066x10^23 - (2.864e17 + 4.225e10 + ... + 0.0

I am trying to use the += operation basically to perform the changes on it (which is depicted in the code below as du[arraySize[1], 1], but it doesn’t seem to actually be affecting that value at all.

function model!(du, u, p, t)
    # Size of the array u.
    arraySize = size(p)

    for i in 1:arraySize[1]
        for j in 2:arraySize[2]
            if p[i, j] != -1
                    
                # Determining the coefficient and value of neighbors of the current node.
                neighbouringValues = getNeighbouringValues(u, i, j)
                neighbouringCoefficients = getInboundCoefficients(p, i, j)
                # We need to calculate the inflow from neighbours to the current node.
                inflow_values = multiply_arrays(neighbouringValues, neighbouringCoefficients)
                inflow = sum(inflow_values)

                # We need to then calculate the outflow from the current node to its neighbors.
                outflow_values = calculateOutflow(p, u, neighbouringCoefficients, i, j)
                outflow = sum(outflow_values)

                # We now sum the changes between the inflow and outflow for the current node.
                du[i, j] = inflow - outflow

                if (i != arraySize[1] && j != 1)
                    # If we are on the right boundary, we need to transfer the monomers from growth to the stockpile.
                    if j == arraySize[2]
                        du[arraySize[1], 1] += p[i, j].Growth * u[i, j]
                        du[i, j] -= p[i, j].Growth * u[i, j]
                    end

                    # We now need to enact changes in Growth/Decay on the stockpile of monomers.
                    # Growing our crystal by size 1 seizes monomers from the stockpile to the crystal.
                    du[arraySize[1], 1] -= outflow_values[1]
                    # Shrinking our crystal by size 1 releases monomers from the crystal to the stockpile.
                    du[arraySize[1], 1] += outflow_values[2]
                end
            end
        end
    end
end

Thank you for your time :slight_smile: !

2

The problem you are running into is that Float64 isn’t anywhere near accurate enough to track this level of accuracy.

3

The issue still occurs whether I use Float64, or the more precise BigFloat type.

4×5 Matrix{BigFloat}:
 -1.0        -1.0            -1.0             -1.0          0.0
 -1.0        -1.0            -1.0              9.89029e-11  1.90176e-17
 -1.0        -1.0          1681.92             0.00290038   1.72268e-10
  6.022e+23   2.86406e+17     4.22592e+10  12482.0          0.000225923
4

I think your logic is wrong. Look at the loops: Everything that happens before i is arraysize[1] does not matter because the very first iteration of j then overwrites du[arraysize[1], 1]. But afterwards you never execute the part the with the += and -= again because the condition

Stays false (because i == arraysize[1]).

Maybe you meant:
if (i != arraySize[1] || j != 1)

EDIT: Just realized that the j loop starts from 2, so the first part of what I wrote is wrong. I was thrown off by the unnecessary check j != 1. Still, you never run the if body again once you reach the row with i==arraySize[1] so there is no subtraction from cells to the right of the lower left one.
There is another issue with your function because you never write values into the first column of du nor any cell that corresponds to a -1. I think that your are required to set every value of du because there are no guarantees what values are present when your function is called by DifferentialEquations. Maybe include a fill!(du, 0) at the top to ensure zeroing it out.