A PDE equation giving meaningless results

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;

Answer 1

[Torsten]

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;

[Dursman Mchabe]

It worked, thanks a lot Torsten. I appreciate.