Parabolic solver of PDE and coefficient in function form.

I want to solve non-stationary parabolic equation d∂u∂t−∇·(c∇u)+au=f , where a=f(u) , but Matlab give out such errors:

Undefined function or variable 'u'.

Error in MoePDE (line 17)
u=parabolic(u0,tlist,b,p,e,t,c,acoeff(p,t,u,tlist),f,d);

where

function coeff = acoeff(p,t,u,time)
% Triangle point indices
it1 = t(1,:);
it2 = t(2,:);
it3 = t(3,:);
% Find centroids of triangles
xpts = (p(1,it1) + p(1,it2) + p(1,it3))/3;
ypts = (p(2,it1) + p(2,it2) + p(2,it3))/3;
uintrp = pdeintrp(p,t,u); % Interpolated values at centroids
A1=-1.0174E-10; A2=2.10708E-10; x0=423976; dx=873999.05251;
Eps1=6.8*8.85e-12;Eps2=22.6*8.85e-12;h=0.002;
coeff=appr2(uintrp,A2,A1,x0,dx,h)/Eps1;
end

and

function v=appr2(u,A2,A1,x0,dx,h)
v = A2 + (A1-A2)./(1 + exp((u-x0)/dx));
end

full code:

vs=1e-12;
c=(-vs*h)/Eps1;
f=0;
d=1;
[p, e, t] = initmesh(g, 'Hmax', 0.08);
pdemesh(p,e,t); axis equal
n=200; 
u0=1600;
tlist=linspace(0,5,260);
u=parabolic(u0,tlist,b,p,e,t,c,acoeff(p,t,u,tlist),f,d);
pdesurf(p,t,u); grid on;

boundary conditions are u=0 on each edge and geometry model(square 1x1) i have export from PDE TOOL.

If I use

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

then result of 'a' will be constant, am i right? (i just build two plots: 1st with a=0.29 and 2nd like@acoeff, and they were same), but i want that a is variable coefficient which varies during the execution of the solver.

Function coefficient 'a' is calculated as written in Documentation:

Calculate Coefficients in Function FormX- and Y-Values

The x- and y-values of the centroid of a triangle t are the mean values of the entries of the points p in t. To get row vectors xpts and ypts containing the mean values:

% Triangle point indices
it1 = t(1,:);
it2 = t(2,:);
it3 = t(3,:);
% Find centroids of triangles
xpts = (p(1,it1) + p(1,it2) + p(1,it3))/3;
ypts = (p(2,it1) + p(2,it2) + p(2,it3))/3;

Interpolated u

The pdeintrp function linearly interpolates the values of u at the centroids of t, based on the values at the points p.

uintrp = pdeintrp(p,t,u); % Interpolated values at centroids

The output uintrp is a row vector with the same number of columns as t. Use uintrp as the solution value in your coefficient calculations.

____________________________________

Thanks in advance for your help

[Torsten]

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

function a = acoeff(p,t,u,time)
a = ...;
end

[Pavel M]

Torsten, thank you for your answer, but I am not up to the end sure that using @acoeff i will get variable coefficient 'a' . In my case a=f(u) and i need to get value of 'a' at different values of 'u'. Maybe, i badly understand how solver 'parabolic' works. Does it call function @acoeff every time when next value of 'u' is calculated?

Thank you in advance

[Torsten]

You can define the "a"-coefficient in function "acoeff" as it depends on u.

[Pavel M]

clc,clear
Eps1=6.8*8.85e-12;Eps2=22.6*8.85e-12;h=0.002;vs=1e-12;
v=1.6e-11/Eps1
b=[1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;48 48 48 48;48 48 48 48;49 49 49 49;48 48 48 48];
g=[2 2 2 2;-0.3 0.3 0.3 -0.3;0.3 0.3 -0.3 -0.3;0.3 0.3 -0.3 -0.3;0.3 -0.3 -0.3 0.3;0 0 0 0;1 1 1 1];
[p, e, t] = initmesh(g, 'Hmax', 0.075);
pdemesh(p,e,t); axis equal
[p1, e1, t1] = refinemesh(g,p,e,t);
% pdemesh(p1,e1,t1); axis equal
x=p(1,:)';
y=p(2,:)';
u0=6600;
n=260;  
tlist=linspace(0,10,260);
a=@acoeff;
% a=0.29
d=1;
f=0;
c=(-vs*h)/Eps1;
u=parabolic(u0,tlist,b,p,e,t,c,a,f,d);

function

function coeff = acoeff(p,t,u,time)
uintrp = pdeintrp(p,t,u); % Interpolated values at centroids

A1=-1.0174E-10;
A2=2.10708E-10;
x0=423976;
dx=873999.05251;

Eps1=6.8*8.85e-12;
% Eps2=22.6*8.85e-12;
h=0.002;

coeff=appr2(uintrp,A2,A1,x0,dx,h)./Eps1;

end

function

function v=appr2(U,A2,A1,x0,dx,h);
v = A2+(A1-A2)./(1+exp.((U-x0)./dx));
end

matlab fives out errors:

Error using  / 
Matrix dimensions must agree.

Error in appr2 (line 2)
v = A2 + (A1-A2)/(1 + exp((U-x0)/dx));

Error in acoeff (line 13)
coeff=appr2(uintrp,A2,A1,x0,dx,h)./Eps1;

[Torsten]

Error message and code are not consistent.

I can't find the line

v = A2 + (A1-A2)/(1 + exp((U-x0)/dx));

in your code, but only

v = A2+(A1-A2)./(1+exp.((U-x0)./dx));

In the last line, you will have to remove the "." after exp

[Pavel M]

Yes, i dont know why error message and code are not consistent.

but v = A2+(A1-A2)./(1+exp((U-x0)./dx)) also dont work

[Torsten]

What do you mean by "don't work" ?

[Pavel M]

That function begin work, when i start other programm

A1=-1.0174E-10; A2=2.10708E-10; x0=423976; dx=873999.05251; Eps=6.3720e-11; h=0.002;
Eps1=6.8*8.85e-12;Eps2=22.6*8.85e-12;vs=1e-12;
v=1.6e-11/Eps1
[T,U] = ode45(@(T,U) (-U*appr2(U,A2,A1,x0,dx,h))/Eps,[0.1:1:100],2800)

After this it can start in u=parabolic...

At the same time, it

function v=appr2(U,A2,A1,x0,dx,h);
v = A2+(A1-A2)./(1+exp((U-x0)/dx));
end

[Torsten]

?

[Pavel M]

when at first i run programm MoePDE, matlab gives error:

Error using  / 
Matrix dimensions must agree.
Error in appr2 (line 2)
v = A2 + (A1-A2)/(1 + exp((U-x0)/dx));

where

 v = A2+(A1-A2)./(1+exp((U-x0)/dx)).

But, then i run other programm, where also i use this function, for example:

A1=-1.0174E-10; A2=2.10708E-10; x0=423976; dx=873999.05251; Eps=6.3720e-11; h=0.002;
Eps1=6.8*8.85e-12;Eps2=22.6*8.85e-12;vs=1e-12;
v=1.6e-11/Eps1
[T,U] = ode45(@(T,U) (-U*appr2(U,A2,A1,x0,dx,h))/Eps,[0.1:1:100],2800)

and after this programm MoePDE also begin work though in function code nothing exchanged

function v=appr2(U,A2,A1,x0,dx,h);
v = A2+(A1-A2)./(1+exp((U-x0)/dx));
end

[Torsten]

Check which function files MATLAB uses on execution.

It can't happen that MATLAB uses

v = A2+(A1-A2)./(1+exp((U-x0)/dx))

to calculate v if it states that there is an error in the line

v = A2 + (A1-A2)/(1 + exp((U-x0)/dx))

(except for the case that this line appears more than once in your code).

[Pavel M]

Maybe it happen because of i use ''function in function''?

function coeff = acoeff(p,t,u,time)
uintrp = pdeintrp(p,t,u); % Interpolated values at centroids
A1=-1.0174E-10;
A2=2.10708E-10;
x0=423976;
dx=873999.05251;
Eps1=6.8*8.85e-12;
% Eps2=22.6*8.85e-12;
h=0.002;
coeff=appr2(uintrp,A2,A1,x0,dx,h)./Eps1;
end