pdepe help! (boundry condition in PDE solver)

40 次查看(过去 30 天)
Hello, I need some help with differential equation solving in MATLAB.
I need to slove heat-conduction equation: Cρ*dT(x,t)/dt = df(dT(x,t)/dx)/dx.
Boundry condition: T(x,0)=T0; dT(L,t)/dx=0; T(0,t) = Tg(t) [is the upper boundary condition and, (here, Tg is an instrument-recorded temperature)]
However, the boundary conditions Tg(t) are not periodic(Temperature is changed with time and have no fixed transform relation with time).
How can use pdepe to solve it?
Can someone help me??
This is my code:
  1 个评论
Yizhou Du
Yizhou Du 2019-1-16
This is my code;
function pdex1
m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
C = 4.0;
k = 419;
density = 8960;
TT = [25 30 35 40 45 50 55 60 65];
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
u1 = sol(:,:,1);
figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = pdex01pde(x,t,u,DuDx)
c = density*C;
f = k * DuDx;
y = u(1);
F = 0*y;
s = F;
% % --------------------------------------------------------------
function u0 = pdex1ic(x);
u0 = 25;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul-TT(i);
ql = 0;
pr = 0;
qr = 1;

请先登录,再进行评论。

采纳的回答

Josh Meyer
Josh Meyer 2019-2-15
Your boundary condition function includes the line
pl = ul-TT(i);
TT(i) is not defined within view of the function, so that is what is causing the error. To make that work you need to pass TT in as an extra parameter. But, still, the index i is undefined. So I would recommend just doing something like:
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
Tg = 25; % Or Tg = h(t) if you can express it as a function of t
pl = ul-Tg;
ql = 0;
pr = 0;
qr = 1;
If you really want to pass in TT as a parameter, change the header of the boundary function to
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t,TT)
and then define the variable
bc = @(xl,ul,xr,ur,t) pdex1bc(xl,ul,xr,ur,t,TT);
This allows pdepe to pass in the 4 inputs it expects, and allows you to pass in the extra TT input from the workspace. You'll need to define the index i somewhere within pdex1bc as well. The call to pdepe then becomes
sol = pdepe(m,@pdex1pde,@pdex1ic,bc,x,t);
Also, although these don't cause errors some other possible issues are:
  • The constant k appears in your flux term f = ..., but not in your original equation.
  • Your source term should just be s = 0 from the original equation, but you have lines to define y, F, and s.
  3 个评论
Cheng Pan Cheng
Cheng Pan Cheng 2019-11-19
Hi, I am a new learner of matlab. I am wondering the meaning of pl, pr, ql, qr, ul, etc. Would you mind explaining it to me? Thanks a lot.
Josh Meyer
Josh Meyer 2019-11-19
The boundary conditions are given as coefficients in an equation that is satisfied on the boundary. The equation is:
To implement the boundary conditions for pdepe, you need to write a function that gives the coefficients p and q in this equation. So, on the left boundary, the boundary equation has certain values for p and q, which are called and . Likewise on the right boundary you have and :
Right:
Left:
The boundary function uses the form
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
In other words, you have access to the values for x, u, and t on each boundary, and need to give the values for p and q on each boundary in terms of those variables.
As an example, let's say on the left boundary you have the condition . In terms of the boundary equation, you have and , since substituting these coefficient values into the boundary equation gives . Now, let's say on the right boundary you have . Since there is still no flux term, you again have , however the expression for p becomes , since if you put these coefficients into the boundary equation you get
One thing to note is that the flux term for f is the same one that is defined in your main equation function for pdepe. It typically contains a partial derivative, but can also have other terms. So generally you'll have q=0 for any boundary that doesn't have a partial derivative in the condition. If your flux term was , but your left boundary condition is , then your coefficients would be and in order to cancel out the 1/2 in the flux term.

请先登录,再进行评论。

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by