Laplace Equation with pdepe command

Question

Thanasis Hou 2017-3-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/330647-laplace-equation-with-pdepe-command

评论： MEENAL SINGHAL 2020-6-14

Hi, I am learning to solve PDEs with pdepe command. First, I try to solve a Laplace equation (1D) setting zero the right coefficients in the PDE. I am expecting the solution to be e.g. y=m*x + b. However I get curves and not lines and I do not know why. I attach the M-files below.

Thanks.

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

请先登录，再进行评论。

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

Answer 1

Torsten 2017-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/330647-laplace-equation-with-pdepe-command#answer_259474

Setting s=1, the solution is of the form u(x)=-1/2*x^2+a*x+b.

Best wishes

Torsten.

6 个评论
显示 4更早的评论隐藏 4更早的评论

MEENAL SINGHAL 2020-6-14

@Torsten

Your suggestion "It's written in the documentation that c identically zero is not allowed. So you could try adding a second (artificial) equation (e.g. du2/dt = 0)." seems valuable to me.

I was wondering what would be the initial condition, if i introduce du2/dt = 0, in addition to my previous equations such that the overall system (a steady state heat transfer equation) remains the same?

Thank you!

-Meenal

MEENAL SINGHAL 2020-6-14

在 MATLAB Online 中打开

pdes.jpg

To add on the complexity, my equations contain discontinuity.

Pdes are defined as in pdes.jpg and the corresponding code is

function forward_compositewalls
global Nci Nco n Nri Nro alphai Betai Gammai alphao Betao Gammao Deltai Deltao et_i et_o 
global Xbw Qi Qo theta_c theta_rad theta_ini Temperature_1 Temperature
Xbw=0.6;
Nci=0; Nco=1;
n=0.5;
Nri=1; Nro=1;
alphai=0.01; alphao=0.01;
Betai=0.04; Betao=0.04;
Gammai=0.05; Gammao=0.03;
Deltai=0.01; Deltao=0.01;
%kratio=0.01;
Qi=1; Qo=1;
theta_c=0.025;
theta_rad=0.025; 
theta_ini=0.03;
et_i=0.6;
et_o=0.4;
L=1;
tend = 3;
m=0;
x = linspace(0,L,21);
t= linspace(0,tend,11);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the solution components as Temperature.
Temperature = sol(:,:,1);
Temperature_1=Temperature';
% A solution profile can also be illuminating.
figure, plot(x,(Temperature_1(:,end)))
% 
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global alphai Betai Gammai alphao Betao Gammao Deltai Deltao Xbw Qi Qo theta_rad 
if x<=Xbw
    c = [0;1];
    f = [1 + Deltai*(u-1).*(DuDx);0];
    s = [Qi.*(1+alphai*u+Betai*u.^2+Gammai*u.^3);0];
else
    c = [0;1];
    f = [1 + Deltao*(u-theta_rad).*(DuDx);0];
    s = [Qo*(1+alphao*u+Betao*u.^2+Gammao*u.^3);0];
end
function u0 = pdex1ic(x)
u0 =[ 0; 0 ];  % i don't know this (What should be the value for the steady case?)
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global Nci Nco n Nri Nro et_i et_o  theta_c theta_rad theta_ini 
pl=[Nci*((ul-theta_c)/(theta_ini-theta_c))^n*(ul-1) + Nri*(1+et_i*(ul-1))*(ul^4-1^4);0];
ql =[-1;0];
pr =[Nco*((ur-theta_c)/(theta_ini-theta_c))^n*(ur-theta_c) + Nro*(1+et_o*(ur-theta_rad))*(ur^4-theta_rad^4);0];
qr=[1;0];