Is it possible to use a variable coefficient "a" in PDE Toolbox with parabolic solver?

Hi all,

I will try to explain my doubt in the better way possible:

I'm trying to solve the reaction-diffusion equation with PDE Toolbox (Matlab) with non-constant coefficients, the syntax to get the solution (u) is:

parabolic - Solve parabolic PDE problem

This MATLAB function produces the solution to the FEM formulation of the scalar PDE problem: 
u1 = parabolic(u0,tlist,b,p,e,t,c,a,f,d)|

c,a,f,d are the coefficients of the parabolic equation:

    d(∂u/∂t)−∇⋅(c∇u)+au=f,

But I want something like this:

u=parabolic(u0,tlist,b,p,e,t,c,@coeffunction,f,d);

where @coeffunction is a function handle which corresponds to the coefficient "a".

My question is: Anyone knows if is it possible to specify "a" as a variable coefficient which varies during the execution of the solver? i.e. If my problem is defined between two times (t0,tf), I want that "a" takes different discrete values during the process of solve (I mean during the execution of parabolic).

Thanks

Answer 1

[Alan Weiss]

The answer is, it depends on your MATLAB version. According to the release notes, in R2012b and later, the a coefficient can me a function of u and its gradient, as well as of x and y. The documentation describes how to represent your coefficients as string functions, or in functional form.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Pablo Díaz]

Hi Alan, thanks for your fast reply!

I'm using Matlab R2015a, yesterday I was able to do this by writing explicitly the coefficients with help of a if-elseif conditional, doing something like this:

if time < time_vector(1)
   a = 'here the value in this time'
elseif time < time_vector(2)
   a = 'here the value in time between time_vector(1) and time_vector(2)'
elseif ...
end

but I have problems when I try to do the same with the help of a function that read "a" during the execution, i.e. I build "a" as a cell array, where the components are the values at the desired times, and when I try to read this array with the help of a for loop during the execution, the program fails "a is not assigned during the execution" ... :(

Btw, do you knows if the solution calculated by the parabolic solver depends on the number of elements which have "tlist"? When I put "a" as a coefficient which depends on the value of u at different times, i.e as a cell array "a" which have the following components:

   a(1) = 'constants.*u(:,time1)'
   a(2) = 'constants.*u(:,time1) + constants.*u(:,time2)'
   ...
   a(n) = 'constants.*u(:,time1) + .... + constants.*u(:,time(n))'

I'm facing the following problem: if I put more elements inside "tlist" but with the same limits (t0 and tf), the solution is not the same, I found a not negligible difference between both solutions... but my knowledge says that a solution mustn't be dependent on the number of elements that have "tlist" while the limits are the same.

Thanks a lot for your helpful help

[Alan Weiss]

I think that you have a few misapprehensions:

• tlist does NOT contain the only times that parabolic might sample. It contains the range of times that are sampled, and the solution times that are reported, but the solver might sample any point in between the first and last points in tlist. So your coefficient function has to be defined for any time in that range.
• a should not be a cell array. You have to give a in the form that is documented, meaning a scalar for a single PDE, or the correct form in case you have a system of PDEs.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Pablo Díaz]

Once again, thank you very much for your help.

Probably, because my bad level of written english, I had explained badly my problem.

• Yes, I know that tlist is only a array that I can use to sample the solution, I know that the solver runs along the "continuous" time (I would like to know if there is any time discretization inside the solver. Computation must be discrete, not continuous)
• Here I could ilustrate what I have:

function a = coeffunction(p,t,u,time)
    global tempo tlist n 
    [a1] = generate_sf(n,tempo, tlist);
    if time <= tempo(1)
        a = '15./(u+2.5)';
        if abs(time-0.1) < 1e-2
            time
            a
            disp('im here tempo 1!')
        end
    elseif time <= tempo(2)
        a = char(a1{1});
        if abs(time-(tempo(1)+10000)) < 10000
            time
            a
            disp('im here tempo 2!')       
        end
    elseif time <= tempo(3)
        a = char(a1{2});
        if abs(time-(tempo(2)+10000)) < 10000
            time
            a
            disp('im here tempo 3!')
        end  
    elseif time <= tempo(4)
        a = char(a1{3});
        if abs(time-(tempo(3)+10000)) < 10000
            time
            a
            disp('im here tempo 4!')
        end    
    elseif time <= tempo(5)
        a = char(a1{4});
        if abs(time-(tempo(4)+10000)) < 10000
            time
            a
            disp('im here tempo 5!')
        end        
    else
        a = char(a1{5});
        if abs(time-(tempo(5)+10000)) < 10000
            time
            a
            disp('im here tempo 6!')
        end
    end

As you can see, to represent the coefficient "a" I have a function "coeffunction" which have a dependency on time and u, I tried to follow the rules of the manual to do this.

Inside "coeffunction" I calculate a cell "a1" with the help of the auxiliar function "generate_sf", finally I have "a1" which has 5 components:

>> a1

a1 = 

      {1x1 cell}    {1x1 cell}    {1x1 cell}    {1x1 cell}    {1x1 cell}

each one component can be converted to a string using "char(a1{i})"

>> a1{1}

ans = 

      '15./(u+2.5).*(1*(exp(-((0.2*2*(u(:,1)*2.5+3.3))./(2.5*(u(:,1)+3.3))+(0.02*2*...'

So, finally I have a "a" coefficient stored in a "cell array" which I access and convert to a "string" element by element. It is wrong?

But the problem that I facing is: Although the times t0 (initial time) and tf (final time) I use to generate the vector "tlist" are the same, If I change the number of elements (for example, tlist1 has 2000 elements, and tlist2 has 4000 elements, buth both had the same initial and final time) the solution is not the same!Even, there are a great discordance between solutions! Why "u" depends on the number of elements in this case? What am I doing wrong?

Thanks a lot for your valuable help

Answer 2

[Pablo Díaz]

In this image, I show the results... this curves corresponds to the same solution of the parabolic PDE problem, but the red line has more elements inside tlist than the red line. Why is this happening?