Laplace Equation with pdepe command

5 次查看（过去 30 天）
Thanasis Hou
Thanasis Hou 2017-3-19
评论： 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.

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

采纳的回答

Torsten
Torsten 2017-3-20
Setting s=1, the solution is of the form u(x)=-1/2*x^2+a*x+b.
Best wishes
Torsten.
  6 个评论
显示 4更早的评论
MEENAL SINGHAL
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
MEENAL SINGHAL 2020-6-14
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];
I am stuck with the initial boundary and need help on this. Moreover, the discontinuity part is creating problem. Any help would be appreciated.
Thanks!
-Meenal

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by