Laplace Equation with pdepe command
显示 更早的评论
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
2017-3-20
0 个投票
Setting s=1, the solution is of the form u(x)=-1/2*x^2+a*x+b.
Best wishes
Torsten.
6 个评论
Thanasis Hou
2017-3-20
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).
From the documentation:
In Equation 1-3, f (x,t,u,∂u/∂x) is a flux term and s (x,t,u,∂u/∂x) is a source term. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c (x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if those values of x are mesh points. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.
Best wishes
Torsten.
Torsten
2017-3-21
If your equation is simply elliptic, use "bvp4c" instead of "pdepe".
Best wishes
Torsten.
Thanasis Hou
2017-3-23
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
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 个）
类别
在 帮助中心 和 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 
