Specify time-dependent PDE coefficients

Hello,

I would ask help concerning the analysis of a reaction-diffusion system with Matlab, with two coupled PDEs. In my problem, I have a defined function for the temperature T(t) (mixed ramps and constants values) that I don't need to solve just because it is imposed. In the PDE system that I need to solve there are two parameters k1 and k2 (entering the coefficient a) that are dependent on T (and then on t, because of T is imposed and not solved as PDE). I have tried to define two functions and then to add a as vector a = [k1(t); k2(t)] as follow:

k1 = @k1f;         % First external function
k2 = @k2f;         % Second external function
m = 0;             % Second derivative null                 
d = [1;1];         % First derivative 1               
c = [0.2;0.1];     % Diffusion              
a = [k1(t);k2(t)]; % Reaction term    
f = [0;0];         % No generation                 
pdec = specifyCoefficients(pdem, 'm', m, 'd', d, 'c', c, 'a', a, 'f', f);

but it doesn't work.

In the PDETool it is simple to do this because I have just to add the function in the user interface and it generates the command pdeseteq(TYPE,C,A,F,D,TLIST,U0,UT0,RANGE).

I would like to know if I can define this parameter without using the PDETool.

Thank you

Answer 1

[Alan Weiss]

I am not sure what you are trying to do. Is the "t" that you "impose" equal to the time in your PDE? Or is there an outer loop that we cannot see that sets "t" for some other purpose?

Also, while you say that you are setting the c coefficient to be a function of "t", it seems in your code that you are setting a as a function of "t". I am going to assume that you made a mistake in setting a, and that you really meant to set c.

If "t" is the time in your PDE, then you have some syntax errors. For a c coefficient that is a function, you need to define the function as

function c = ccoeff(region,state)

It seems that you want a "short" 2-element c output. But, as the documentation clearly states, the returned c is a matrix of size N1-by-Nr, where N1 = 2 for your case, and Nr is the number of elements in each region structure that MATLAB passes to your solver. So your c coefficient function should be something like this:

function c = ccoeff(region,state)
t = state.time;
k1 = k1f(t);
k2 = k2f(t);
c = repmat([k1;k2],1,length(region.x)); % 2-by-Nr matrix c

When you set coefficients, set c as @ccoeff.

I hope that this enables you to continue with your work.

Alan Weiss

MATLAB mathematical toolbox documentation

[MG]

Dear Alan,

thank you for your kind reply.

At first, to reply to your questions: - Is the "t" that you "impose" equal to the time in your PDE? Or is there an outer loop that we cannot see that sets "t" for some other purpose? It is the same in my PDEs.

- Also, while you say that you are setting the c coefficient to be a function of "t", it seems in your code that you are setting a as a function of "t". I am going to assume that you made a mistake in setting a, and that you really meant to set c. No, I meant a instead of c, sorry (I will edit the main post).

I have tried to modify my code following your suggestions but, unfortunately, it doesn't work. By running the new code I obtain the same error:

Function handle specifying a coefficient must accept two
input arguments and return one output argument.

Also for 'a' coefficient I need a N1-by-Nr matrix and I have verified that the repmat command returns the correct result.

I don't know what is the problem but I will try to work on this new code and understand more in deep.

If you have other suggestions, these are warmly welcome.

Thank you.

[Alan Weiss]

While you didn't show your function for specifying the a coefficient, the error is clear: the function needs to accept region and state inputs, and return a single matrix output. I suggest that you ensure that your function does exactly that.

For more help, if necessary, please show BOTH the coefficient function AND the line of code where you call specifyCoefficients.

Alan Weiss

MATLAB mathematical toolbox documentation

[MG]

This is the temperature that I would impose:

function Tc = Tcf(~,state)

t = state.time;

Tc = (ramp(t,170,0) + ramp(t,-170,0.5) + ramp(t,-170,1.5) + ...
    ramp(t,170,2) + ramp(t,170,3.5))';

end

while k1 is defined as:

function k1 = k1f(~,state)
E = 130e3;  % J/mol
A = 3.8e3;   % 1/h
R = 8.314;   % J/mol K

t = state.time;
Tc = Tcf(t);

k1 = A*exp(-E./(R*Tc));

end

The same is for k2, which has different E and A. The 'a' coefficient is:

function a = acoef(region,state)

nr = numel(region.x);
t = state.time;

k1 = k1f(t);
k2 = k2f(t);

a = repmat([k1;k2],1,nr);

end

and finally I use this command:

pdec = specifyCoefficients(pdem, 'm', m, 'd', d, 'c', c, 'a', @acoef, 'f', f);

where m, d and c are those previously defined.

I have finished my ideas to fix it...

[Torsten]

In Tc, k1 and k2:

t=state

instead of

t=state.time

Best wishes

Torsten.

[MG]

Same error...

Function handle specifying a coefficient must accept two
input arguments and return one output argument.

What a mess...

[Alan Weiss]

Thank you for supplying the details.

Functions Tc and k1f are NOT functions of two input variables. They are functions of ONE variable, usually called t. Each should be written as something like

function k1 = k1f(t)
% do NOT put a call to region or state here

Similarly,

function Tc = Tcf(t)
% These are functions of one input only

One more thing: did you write the ramp function? Or are you calling the Simulink ramp block? I am asking because that part of your function is something that I still do not understand. I see that you have a transpose as the last operation in the function when you call it in Tc, but I would expect the function to return a scalar value, so I am suspicious. If you are calling Simulink, well, that is yet another problem.

Your acoef function looks good to me.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[MG]

Thank you Alan,

I corrected Tc, k1 and k2 as you suggested. Still doesn't work. However, I have this warning now:

Warning: Failure at t=1.000000e+00.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (1.776357e-15) at time t.

I should be able to manage this that is better than the previous.

Concerning the ramp function, I have used a pre-compiled script from the File exchange section. It is the follow:

function r = ramp(t, slope, startTime)

N = numel(t);

r = zeros(N,1);

if median(diff(t))>0
    startInd = min(find(t>=startTime));
    popInd = startInd:N;
elseif median(diff(t))<0
    startTime = - startTime;
    startInd = max(find(t>=startTime));
    popInd = 1:startInd;
    slope = -slope;
end

r(popInd) = slope.*(t(popInd) + startTime) - 2*startTime*slope;

return

I don't like it so much for this particular application, I tried to use a if cycle like this:

function Tc = Tramp(t)

if t < 0.5
    Tc = 170*t;
elseif t < 1
    Tc = 1;
elseif t < 1.5
    Tc = -170*t;
elseif t < 2
    Tc = -1;
else
    Tc = 170*t;
end

But it returns only the part of the function defined in the else. I don't know why.

I could try to change the Tcf function to fix the script.

So, thank you Alan and Torsten for your useful suggestions.

Bests