solving a 200 by 200 by 200 array ode

Hi

I am trying to model the heat transfer of a 3D object obtained by ct scan, which is a 200 by 200 by 200 array. when I use ode solver to solve that, I am out of memory. is there any way I can bypass this issues ?

thank you very much for your help

[Walter Roberson]

What is your code? What is the sizes of the matrices involved (factoring the array sizes from the image suggests that the array is not 200 x 200 x 200.)

[Huijian Huang]

Hi Walter you are right, it does seen to be something wrong. here is my codes

the array is 112 x 127 x 226 in precise

regards

Answer 1

[Walter Roberson]

You should firmly avoid using global in the ode function. See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

You should vectorize your code. Your code can be theoretically broken into two parts:

1. fully vectorizable case for everything that is completely "inside" the cube
2. 8 boundary cases

As I look at your boundary cases, it looks to me as if you would probably be able to fully vectorize all of your code, if you were to take the original cube of data, and pad it on all sides with a copy of the edge data.. like

                   A A B C D D
   A B C D         A A B C D D
   E * * F  ->     E E * * F F
   G H I J         G G H I J J
                   G G H I J J

If you do that then you can fully vectorize over the entire inside of the augmented volume. This is because the contribution of each term depends upon the difference between adjacent locations, and if you duplicate data like this, the difference between adjacent locations would be 0 at the boundary, leading to the boundary terms contributing 0 to the sums, same as if you had not included them in the calculations.

That said:

The reason that the memory requirements are so high are that gradient estimation must calculate the effect of each variable upon each other variable. With 112 x 127 x 226 locations, that requires calculating (112 x 127 x 226) * (112 x 127 x 226) values. At least two matrices these size are required, so that finite differencing can be performed in order to determine how to change the variables relative to each other.

In your case, only the adjacent variables can affect each other, leading to something that could be represented as a sparse multiband matrix. (Unfortunately because of the 3-space nature of the calculations, you cannot use a simple tridiagonal matrix.)

You might therefore be able to take advantage of one of the methods of specifying sparse jacobian calculations. See https://www.mathworks.com/help/optim/ug/nonlinear-equations-with-jacobian-sparsity-pattern.html that talks about one of the possibilities briefly and goes into more depth on another possibility.

[Huijian Huang]

Thanks Walter

your suggestion and information is very useful. however, I am still an intermediate matlab user, would you be able to kindly show me some example of how to apply your theory. a simple example will be truly appropriated.

best regards

[Walter Roberson]

Unfortunately I have never implemented sparse jacobian or hessian calculation myself.