pdepe boundary condition help!!

H0 = 1;

% Spatial and time mesh

y = linspace(0, H0, 20);

t = linspace(0, 1, 200);

% Solve the PDE

sol = pdepe(0, @pdefun, @icfun, @bcfun, y, t);

% Extract the solution

lambda = sol(:,:,1).';

function [c, f, s] = pdefun(y, t, lambda, DlambdaDx)

nu = 1.7e-28;

X = 0.1;

G0 =1e7;

kB = 1.3806488e-23;

theta = 298;

m_d = 0.5/(kB*nu*theta);

% Define the PDE coefficients

c = 1 / nu;

f = m_d/lambda * (G0*nu*(1+1/lambda^2))*DlambdaDx + kB*theta*(1/(lambda*(lambda-0.99))-1/lambda^2-2*X/lambda^3)*DlambdaDx;

s = 0;

end

function lambda0 = icfun(y)

% Initial condition

if y < 1 %H0

lambda0 = 1;

else

lambda0 = 1.498;

end

end

function [pl, ql, pr, qr] = bcfun(yl, lambda_l, yr, lambda_r, t)

% Boundary conditions

pl = 0;

ql = 1;

pr = lambda_r-1.498 ;

qr = 0;

end

I am trying to solve this problem and confused about the boundary condition where dLambda/dX = 0 at x = 0. I am not getting the results. Can someone explain how to apply initial and boundary condtions? Thank you very much.

Answer 1

[Torsten]

Note that the f you define in "pdefun" doesn't fit the formula from the paper. You forgot to multiply both terms with m_d/lambda, not only the first term.

    H0 = 1;
    % Spatial and time mesh
    y = linspace(0, H0, 40);
    t = linspace(0, 1, 20);
    % Solve the PDE
    sol = pdepe(0, @pdefun, @icfun, @bcfun, y, t);
    % Extract the solution
    lambda = sol(:,:,1);
    plot(y,[lambda(1,:);lambda(5,:);lambda(10,:);lambda(15,:);lambda(end,:)])

    
function [c, f, s] = pdefun(y, t, lambda, DlambdaDx)
    nu = 1.7e-28; 
     X = 0.1;
    G0 =1e7; 
    kB = 1.3806488e-23; 
    theta = 298;
    m_d = 0.5/(kB*nu*theta);
    
    % Define the PDE coefficients
    c = 1 / nu;
    f = m_d/lambda * (G0*nu*(1+1/lambda^2) + kB*theta*(1/(lambda*(lambda-0.99))-1/lambda^2-2*X/lambda^3))*DlambdaDx;
    s = 0;
end
function lambda0 = icfun(y)
    % Initial condition
    if y < 1 %H0
        lambda0 = 1;
    else
        lambda0 = 1.498;
    end
end
function [pl, ql, pr, qr] = bcfun(yl, lambda_l, yr, lambda_r, t)
    % Boundary conditions
    pl = 0;
    ql = 1;
    pr = lambda_r-1.498 ;
    qr = 0;
end

[KANCHETI]

Thank you very much for your help. Also could you tell me how can i calculate μ from equation 44 for every value of λ. Also if λ = 1, i need to replace log(1-1/λ) by log(1-0.999/λ).

[Torsten]

Also could you tell me how can i calculate μ from equation 44 for every value of λ.

μ = nu*(G0*λ-G0/λ+ kB*theta/nu*(log(1-1/λ)+1/λ+X/λ^2))

Also if λ = 1, i need to replace log(1-1/λ) by log(1-0.999/λ)

I don't know about the physical background of your problem. So I can't give advice how the case λ = 1 is to be handled in the equations.

[KANCHETI]

Thank you for your time.