Parabolic solver of PDE and coefficient in function form.

1 次查看（过去 30 天）

Pavel M 2019-6-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/465184-parabolic-solver-of-pde-and-coefficient-in-function-form

编辑： Pavel M 2019-6-4

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

12 个评论
显示 10更早的评论隐藏 10更早的评论

Pavel M 2019-6-4

在 MATLAB Online 中打开

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 2019-6-4

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 2019-6-4

编辑：Pavel M 2019-6-4

在 MATLAB Online 中打开

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