PDE toolbox: Undefined function 'mtimes' for input arguments of type 'function_handle'.

Hi everyone,

I have a bit of a hard time finding a way to use the output of my function handle:

effStress = @(~,state) (2700*9.81*4000*1e-6)-state.u;
K = 1e-12-(0.04343*((0.012+0.013)/2));
Dcore=(K*effStress);

I understand that I apparently cannot multiply K with effStress, but even a matrix multiplication doesnt work.

Could anyone help me on that?

Cheers, Flo

[Dennis]

You probably need to provide an input to effStress. Small example:

eff=@(a) a+1; %eff is a function handle
%A=3*eff; %this does not work
A=3*eff(3) %this does.

[Flo brd]

I did try something similar. It didn't work. My issue is that I am working with the PDE toolbox. I am trying to apply a non-constant variable (with the C coefficient if that helps)

Answer 1

[Stephen23]

A function handle is not a numeric array, it cannot be multiplied. You could either

1. evaluate the function to get a numeric value, and multiply that value:

>> S.u = 4;
>> effStress = @(~,state) (2700*9.81*4000*1e-6)-state.u;
>> K = 1e-12-(0.04343*((0.012+0.013)/2));
>> K*effStress(0,S) % evaluate effStress to get an output.
ans = -0.055345

2. define a new function that calls your function:

>> fun = @(s)K*effStress(0,s); % does not evaluate effStress.
>> fun(S)
ans = -0.055345

[Flo brd]

Thanks for your answer. My problem is that I don't really have an idea of what happens in effStress, as I am using the PDE toolbox. I'd just like to get the state of my function u through this function handle. Other example of badly assessed handles have been solved like this one:

https://uk.mathworks.com/matlabcentral/answers/345324-pde-toolbox-how-to-specify-coefficient-for-temperature-dependent-properties

But none of that seems to work in my case...

[Stephen23]

"My problem is that I don't really have an idea of what happens in effStress, as I am using the PDE toolbox"

It is not clear what you mean: a function handle is a function handle, regardless of what toolboxes you might have installed. You can evaluate them, pass them as variables, or use them to define other function handles. What do you want to do with your function handle?

"I'd just like to get the state of my function u through this function handle"

Then you might want to evaluate the function handle to get an output. But it is not clear, because the "state" of your function is a property of your data, not a property of any function handle (i.e. the state is determined by some meaning you attribute to some data. How this "state" relates to any function handle depends on your algorithm and how you have written your code... not things that I know).

[Ravi Kumar]

Can you provide the context of usage with a complete example?? Which function of PDE Toolbox are you referring to? Is your question, what is the value of state.u that solver passes when solve a model?

[Flo brd]

Sorry for the lack of data I was able to provide in the first place. So as said earlier, I am working with the PDE toolbox. I am trying solve for the equation of fick, and that works well I guess. Now I need to update a non-constant variable as the PDE solver solves the partial differential equation.

Here is the code we are interested in

N=1;
model = createpde(N);

%define boxes et generate mesh
height = 10;
length = 30;

r1 = [3 4 0 10 10 0 0 0 height height];
r2 = [3 4 10 20 20 10 0 0 height height];
r3 = [3 4 20 length length 20 0 0 height height];
gdm = [r1; r2;r3]';
g = decsg(gdm,'R1+R2+R3',['R1'; 'R2'; 'R3']');
geometryFromEdges(model,g);
msh =generateMesh(model,'Hmax',2);
nodes = msh.Nodes;

%set id to nodes
nodes(3,:)= linspace(1,size(nodes,2),size(nodes,2));

% BC
applyBoundaryCondition(model,'dirichlet','edge',[7,8,9,4,5,6],'h',0);
applyBoundaryCondition(model,'dirichlet','edge',[3,10,1,2],'h',1);

%specify coefficients for the equation described in the PDE
a =0;
f = 0;
d = 1;

Dcore= @(location, state) Diffusivity(location, state, 1e-16,6000, 2700, 200);
specifyCoefficients(model,'m',0,'d',d,'c',Dcore,'a',a,'f',f, 'face',1);

Followed by the Diffusivity function:

function Dh = Diffusivity(location,state,ki,depth,density,T)

      N = 1;                  % Number of equations
      nr = length(location.x);  % Number of columns
      Eff = zeros(N,nr);     % Allocate to avoid size issue

      confiningStress = density*9.81*depth;
      effStress= (confiningStress*1e-6) - state.u;

      tmp = effStress .*Eff;

      Cp =  0.000000001;
      T = T+273.15;
      A=2.414*10^-5;
      B=247.8;
      C=140;

      %temperature dependent viscosity
      u = A*10^(B/(T-C));

      %pressure dependent k
      gamma = 0.1;
      K=ki-(0.2*gamma.*tmp);

      Dh = K./(u*Cp);
  end

that compiles, but I am really unsure I've understood the way the c coefficient works https://uk.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html#bu5xabm-1

In my understanding, state.u represents the state of the equation u at a certain point of the solving process. As my function depends on a value u, I am trying to get that value through state.u to update Dcore accordingly.

Am I right in my understanding of function handle?

I've tried using those threads as a base reflection:

https://uk.mathworks.com/matlabcentral/answers/383509-non-constant-parameters-for-pde-thermal-solver

https://uk.mathworks.com/matlabcentral/answers/345324-pde-toolbox-how-to-specify-coefficient-for-temperature-dependent-properties

Cheers, Flo

[Flo brd]

I changed a bit my function, following what I've understood from the matlab help (and I am pretty my understanding is wrong), and my error changed to a warning. it doesn't compile but maybe you help me on that one too. here is my function:

function Dh = Diffusivity(location,state,ki, depth,density, T)

      N = 1;                  % Number of equations
      nr = numel(location.x);  % Number of columns
      Dh = zeros(N,nr);     % Allocate to avoid size issue
  %     Eff(1,:) = ones(1,nr);

      confiningStress = density*9.81*depth*1e-6;

      matConf = zeros(1,nr); 
      matConf(1,:)=confiningStress;

      Cp =  0.000000001;

      T = T+273.15;
      A=2.414*10^-5;
      B=247.8;
      C=140;

      %temperature dependent viscosity
      u = A*10^(B/(T-C));

      %pressure dependent k
      gamma =0.1

     Dh(1,:)=(ki-(0.2*gamma.*(matConf(1,:) - state.u(1,:))))./(u*Cp);

end

the error is now

Warning: Failure at t=1.132789e-12. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.231174e-27) at time t.

I am completely lost.

Flo