A PDE equation giving meaningless results

조회 수: 1 (최근 30일)
Dursman Mchabe
Dursman Mchabe 2018년 9월 19일
댓글: Dursman Mchabe 2018년 9월 19일
Hi Everyone,
I am trying to solve a simpe PDE equation using pdepe, however I get meaningless results. Please see the code below. I don't know where could I be going wrong. The equations are:
% simple PDE % Equation: dC/dt = K * d^2C/dx^2 - kr * C
% IC: C0 = 0.762471166
% BC: x = 0, t > 0, C = Cfinal
% BC: x = delta, t > 0, C = 0
The code looks like this:
function SimplePDE
% simple PDE
% Equation: dC/dt = K * d^2C/dx^2 - kr * C
% IC: C0 = 0.762471166
% BC: x = 0, t > 0, C = Cfinal
% BC: x = delta, t > 0, C = 0
global K kr kl dp delta
K = 1.1e-8;
kr = 4.16e-5;
kl = 9.93e-3;
dp = 3.64e-5;
delta = (dp * K)/(2 - kl * dp);
m = 0;
x = linspace(0,delta,100);
t = linspace(0,30,100);
sol = pdepe(m,@DRpde,@DRic,@DRbc,x,t);
C = sol(:,:,1);
figure, plot(t,C(:,1),'r--')
xlabel('Time (sec)')
ylabel('Conc (mol/m^{3})')
hold on
%%EXPERIMENTAL
time = [1;
3;
5;
7;
9;
11;
13;
15;
17;
19;
21;
23;
25;
27;
29;
];
Conc = [0.762471166
0.630957344
0.10964782
0.051286138
0.027542287
0.020417379
0.014791084
0.010715193
0.009120108
0.007943282
0.00691831
0.006309573
0.005623413
0.005128614
0.004897788
];
plot(time, Conc,'bv', 'markersize',6,'markerfacecolor','b')
xlabel('Time (sec)')
ylabel('Conc H^+ (mol/m^3)')
hold off
% --------------------------------------------------------------
figure, surf(x,t,C)
zlabel('Time (sec)')
xlabel('Conc (mol/m^{3})')
ylabel('Boundary Layer Thickness (m)')
% --------------------------------------------------------------
function [c,f,s] = DRpde(x,t,C,DCDx)
global K kr
c = 1;
f = K * DCDx;
s = -kr * C;
% --------------------------------------------------------------
function C0 = DRic(x)
C0 = 0.762471166;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = DRbc(xl,Cl,xr,Cr,t)
pl = Cl;
ql = 0;
pr = Cr;
qr = 0;

채택된 답변

Torsten
Torsten 2018년 9월 19일
Try
function main
% simple PDE
% Equation: dC/dt = K * d^2C/dx^2 - kr * C
% IC: C0 = 0.762471166
% BC: x = 0, t > 0, C = Cfinal
% BC: x = delta, t > 0, C = 0
global K kr kl dp delta Cfinal
K = 1.1e-8;
kr = 4.16e-5;
kl = 9.93e-3;
dp = 3.64e-5;
delta = (dp * K)/(2 - kl * dp)
Cfinal = 0.3;
m = 0;
x = linspace(0,delta,100);
t = linspace(0,1e-18,8);
sol = pdepe(m,@DRpde,@DRic,@DRbc,x,t);
C = sol(:,:,1);
plot(x,C(1,:),x,C(2,:),x,C(3,:),x,C(4,:),x,C(5,:))
zlabel('Time (sec)')
xlabel('Conc (mol/m^{3})')
ylabel('Boundary Layer Thickness (m)')
% --------------------------------------------------------------
function [c,f,s] = DRpde(x,t,C,DCDx)
global K kr
c = 1;
f = K * DCDx;
s = -kr * C;
% --------------------------------------------------------------
function C0 = DRic(x)
global delta Cfinal
if x==0
C0 = Cfinal;
elseif x==delta
C0 = 0;
else
C0 = 0.762471166;
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = DRbc(xl,Cl,xr,Cr,t)
global Cfinal
pl = Cl - Cfinal;
ql = 0;
pr = Cr;
qr = 0;
  댓글 수: 1
Dursman Mchabe
Dursman Mchabe 2018년 9월 19일
It worked, thanks a lot Torsten. I appreciate.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Eigenvalue Problems에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by