Alternative to 'evalin' for pdepe solver

I am using the built-in matlab 'pdepe' solver to solve a 1D heat transfer problem. I am calling the function like this:

sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

where m is a symmetry variable, x is the spatial discretization, and t is the time discretization. My boundary condition function heatbc and my initial condition function heatic are fixed. However, my pde function heatcyl has material property values that are calculated earlier in the code. The function looks like this:

function [c,f,s] = heatcyl(x,t,u,dudx)

w=evalin('base','w'); % Can't pass these in as arguments when using solver

gamma=evalin('base','gamma');

c = 1;

f = dudx*gamma;

s=funcHj(x)*cos(w*t)*w;

end

For reference, the variable w is the signal frequency, and gamma is the material diffusivity. These change each time the solver pdepe is called. The pde function also calls a function funcHj:

function Hj = funcHj(x)

i_ply=evalin('base','i_ply');

ply_thickness=evalin('base','ply_thickness');

surfthick=evalin('base','surfthick');

if (1+floor((x-surfthick)/ply_thickness)) == i_ply

Hj = 1;

else

Hj = 0;

end

end

Essentially, I have a heating source that changes throughout the thickness of my material, and this Hj function is how I determine the thermal forcing signal at each position along the spatial axis. These will change each time the solver calls on the pde function heatcyl.

As you can see, I am accessing 5 variables (w, gamma, i_ply, ply_thickness, surfthick) within the pde function. Getting these values using Evalin has be working 'fine' for me when I call the pde from the base program. However, I now want to call this solver from within a function. Evalin('base','variable') is not working anymore since the variables are no longer in the 'base' workspace. Additionally, the variables are not in the 'caller' workspace, either. (they are in the workspace of the caller's caller)

For clarity, this set-up works:

%% Start main program

for w = [1,2,3]

sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

end

%% End main program

but this one does not:

%% Start main program

heating = solve_for_frequency(freq1)

%% End main program

function [foo] = solve_for_frequency(w)

sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

foo = sol(:,:,1) % extracting the intersting part of the solution

end

In my search for solutions to my problem, I've found that using Evalin is bad practice, so I'd like to find another way to do this. I've tried saving the variables in a temporary .mat file and loading them from the PDE, but I've read that this is not any better than using Evalin. Additionally, using Load causes the solver to run for 10x longer than using Evalin.

I made an attempt to modify the Mathworks pdepe function to accept my 5 variables and pass them into the pde function, but I was not successful. Additionally, I would like to include my code with a journal publication and I imagine that there will be intelectual property issues preventing me from sharing a modified Mathworks function.

Does anyone have a suggestion on how I can pass my 5 arguments into the pde function?

Many thanks.

Answer 1

[Stephen23]

"Does anyone have a suggestion on how I can pass my 5 arguments into the pde function?"

The recommended, documented, efficient approach is to parameterize your function/s:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

If the parameters were also being changed within the functions you are calling then you would not have much choice but to use nested functions or write a special class for that data. But as your functions merely seem to consume those parameter values, then you can use a simple anonymous function, just as that link above shows.

w = ..;

gamma = ..;

i_ply = ..;

ply_thickness = ..;

surfthick = ..;

fnh = @(x,t,u,d) heatcyl(x,t,u,d,w,gamma,i_ply,ply_thickness,surfthick);

sol = pdepe(m,fnh,@heatic,@heatbc,x,t);

where:

function [c,f,s] = heatcyl(x,t,u,dudx,w,gamma,i_ply,ply_thickness,surfthick)

c = 1;

f = dudx*gamma;

s = funcHj(x,i_ply,ply_thickness,surfthick)*cos(w*t)*w;

end

function Hj = funcHj(x,i_ply,ply_thickness,surfthick)

if (1+floor((x-surfthick)/ply_thickness)) == i_ply

Hj = 1;

else

Hj = 0;

end

end

Note that the content of funcHj can be simply replaced with:

Hj = (1+floor((x-surfthick)/ply_thickness)) == i_ply;

Which hints that most likely funcHj is superfluous anyway:

function [c,f,s] = heatcyl(x,t,u,dudx,w,gamma,i_ply,ply_thickness,surfthick)

c = 1;

f = dudx*gamma;

Hj = (1+floor((x-surfthick)/ply_thickness)) == i_ply;

s = Hj*cos(w*t)*w;

end

[Alex Quinlan]

That's brilliant. That fixed my argument issue, sped up my code, and simplified the Hj function. Thanks a million!