Computation time of PDE Toolbox increases drastically when defining the parameter 'f' as a function handle that uses state.u as input

I use the PDE Toolbox to simulate the diffusion of molecules in a 3D-space (Matlab 2015b). Initially I included separate steps outside the toolbox to account for secretion and degradation of those molecules. However, these processes happen simultaneously, so I tried to implement them into the toolbox using the parameter f.

I created a function ( f = @ffunction ), which calculates the secretion rates and additionally concentration changes due to a predefined half-life after every iteration by accessing state.u. While the first part hardly changed the calculation time, the degradation - using state.u as input - lead to a huge increase. I then found out, that apparently this step alone resulted in the function being called far more often than before (increase from 30 to about 300 times).These are the lines that cause the increase in calculation time:

function f = ffunction(region,state)
ht = evalin('base',ht1);
u = state.u;
tht = 2^(1/ht1);
u1 = u./tht;
f = u1-u;
end

I read a few topics about Matlab integrators requiring smooth functions, so might that be a potential cause?

Answer 1

[Alan Weiss]

I am not sure whether you are using ht1 to represent time, in which case I think that you would do better not to use a global variable ht1, but instead use the state structure state.time. But I could easily be wrong.

When you include a nonlinear function for f, the solver solves a sequence of linear problems, and iterates to attempt to find a fixed point for the nonlinear problem. This is the reason that there are many more function evaluations. For a sketch of the nonlinear solver, see Nonlinear Equations.

Alan Weiss

MATLAB mathematical toolbox documentation

[John Smith]

Thank you for your response. ht1 in this case represents the half-life of the molecule, which I define at the start of my code, so I think I have to use a global variable here.

But if the nonlinearity is the problem here, why does the secretion, which only takes place at distinct coordinates and does not follow a special pattern not cause a similar result?

Apart from that, is there a possibility to obtain an approximation and thereby save some computation time? The concentration is only a contributing factor to my model and not the primary output, so that there is some margin for "error".