Specify time-dependent PDE coefficients

Question

MG 2017-3-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329764-specify-time-dependent-pde-coefficients

评论： MG 2017-3-14

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Weiss 2017-3-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329764-specify-time-dependent-pde-coefficients#answer_258643

在 MATLAB Online 中打开

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

MG 2017-3-14

在 MATLAB Online 中打开

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 2017-3-14

在 MATLAB Online 中打开

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 2017-3-14

在 MATLAB Online 中打开

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