Solving partial differential equation- SOLVED, can be closed

I solved my problem by now, answers are not longer necessary.

Hi all,

I am trying to solve the following partial differential equation:

in which tau,xi and alpha are dimensionless parameters which I can replace by the real parameters(which are shown in the matlab script) My initial condition would be that alpha(x,t=0) = 1 and that

If necessary the boundary condition

d alpha(x,t=0)/dt = 0

can be used too.

I found the function PDEX3 in Matlab and thought that that is useful for my case, but I am struggling with the boundary conditions/initial conditions in this file. I have pasted my code below, I was wondering if anyone knows a more simple way to solve this instead of with this file, or can help me out with the existing? For me the use would be to then apply the real parameters and do some sensitivity analysis with real life values.

Thanks in advance,

Martin

function pdex3
close all;
clear all;
%PDEX3  Example 3 for PDEPE
%   This example illustrates the use of PDEVAL to obtain partial derivatives
%   that are part of solving the problem. A relatively fine spatial mesh is
%   needed to obtain accurate partial derivatives. Values are required at a
%   relatively large number of times. 
%
%   This problem arises in transistor theory. It is used in [1] to illustrate
%   solution of PDEs with series. In the form expected by PDEPE, the single PDE
%   is  
%
%    [1] .*  D_ [u] =  D_ [   d*Du/Dx   ] +  [ -(eta/L)*Du/Dx ]
%            Dt        Dx
%    ---        ---        -------------      ----------------
%     c          u         f(x,t,u,Du/Dx)      s(x,t,u,Du/Dx)
%
%   Here d and eta are physical constants. The equation is to hold on an
%   interval 0 <= x <= L for times t >= 0. For the problem at hand, L = 1. 
%   The initial condition  
%
%      u(x,0) = (K*L/d)*((1 - exp(-eta*(1 - x/L)))/eta)
%
%   involves another physical constant K. The left bc is u(0,t) = 0:
%
%      [u]    +     [0] .* [       Du/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(0,t,u)      q(0,t)        f(0,t,u,Du/Dx)        0
%
%   The right bc is u(L,t) = 0:
%
%      [u]    +     [0] .* [       Du/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(L,t,u)      q(L,t)        f(L,t,u,Du/Dx)        0
%
%   See the nested functions PDEX3PDE, PDEX3IC, and PDEX3BC for the coding of 
%   the problem definition.  
%
%   A quantity of physical interest is the emitter discharge current
%
%      I(t) = (I_p*d/K)*Du(0,t)/Dx
%
%   where I_p is another physical constant. This is meaningful only for t > 0
%   because of the inconsistency in the boundary values at x = 0 for t = 0 and
%   t > 0. This is reflected in the failure to converge of the series solution
%   of [1] for t = 0. Because values for the physical parameters are not
%   specified in [1], nominal values are assumed. I(t) changes rapidly, so
%   output at a relatively large number of t are needed. Recall that the number
%   and placement of the entries of t have little effect on the cost of the
%   computation. Also, the approximation of the solution is only second order
%   in space and the approximation of the partial derivative is of still lower
%   order. Nevertheless, there is reasonable agreement between the computed
%   emitter discharge current and that obtained from the series. 
%
%   [1] E.C. Zachmanoglou and D.L. Thoe,Introduction to Partial Differential 
%       Equations with Applications, Dover, New York, 1986. 
%
%   See also PDEPE, PDEVAL, FUNCTION_HANDLE.

     Lawrence F. Shampine and Jacek Kierzenka
     Copyright 1984-2009 The MathWorks, Inc.
     $Revision: 1.7.4.3 $  $Date: 2009/04/21 03:23:47 $

% Problem parameters, shared with nested functions.
r = 0.03;                       %m
C = sqrt((sqrt(3)-0.5*pi()));
A = (C^2) * (r^2);              %m2
surf_tension = 0.03;            % surface tension, 0.03
rho = 1;                        % kg/m3, water
g = 9.81;                       % gravitational constant
x_0 = sqrt(C*surf_tension/(rho*g));
alpha = A/(x_0^2);              %liquid fraction
viscosity = 0.001;              %liquid viscosity, water
eta_star = 150*viscosity;       
t_0 = eta_star/sqrt(C*surf_tension*rho*g);

tau = t/t0;
L = 1;
d = sqrt(alpha)/2;
eta = alpha;
K = d;
I_p = 1;

m = 0;
x = linspace(0,L,101);
xnew = x./x_0;
t = linspace(0,10,51);
tnew = t./t_0;
sol = pdepe(m,@pdex3pde,@pdex3ic,@pdex3bc,xnew,tnew);
u = sol(:,:,1);

figure;
surf(xnew,tnew,u);
title('Numerical solution computed with 41 mesh points.');
xlabel('Distance x');
ylabel('Time t');

nt = length(t);
I = zeros(1,nt);
seriesI = zeros(1,nt);
iok = 2:nt;
for j = iok
  % At time t(j), compute Du/Dx at x = 0.
  [~,I(j)] = pdeval(m,x,u(j,:),0);
  seriesI(j) = serex3(t(j));
end
% I(t) = (I_p*d/K)*Du(0,t)/Dx
I = (I_p*d/K)*I;

figure;
plot(t(iok),I(iok),t(iok),seriesI(iok));
legend('From PDEPE','From series');
title('Emitter discharge current I(t)');
xlabel('Time t');

    % -----------------------------------------------------------------------
    % Nested functions -- parameters are provided by the outer function.
    %

    function [c,f,s] = pdex3pde(x,t,u,DuDx)
    % Evaluate the differential equations components.
      c = 1;                 
      f = d*DuDx;            
      s = -d*eta*DuDx;       
    end
    % -----------------------------------------------------------------------

    function u0 = pdex3ic(x)
    % Evaluate the initial conditions components.
      u0 = (K*L/d)*(1 - exp(-eta*(1 - x)))/eta;
    end
    % -----------------------------------------------------------------------

    function [pl,ql,pr,qr] = pdex3bc(xl,ul,xr,ur,t)
    % Evaluate the boundary conditions components.
      pl = ul;                 
      ql = 0;                 
      pr = ur;                 
      qr = 0;                 
    end 
    % -----------------------------------------------------------------------

    function It = serex3(t)
    % Approximate I(t) by 40 terms of a series expansion.
      It = 0;
      for n = 1:40
        temp = (n*pi)^2 + 0.25*eta^2;
        It = It + ((n*pi)^2 / temp)* exp(-(d/L^2)*temp*t);
      end
      It = 2*I_p*((1 - exp(-eta))/eta)*It;
    end
    % -----------------------------------------------------------------------  

end  % pdex3

[Bill Greene]

If you can state your boundary conditions, someone might be able to help you. You say,

d alpha(x,t=0)/dt = 0

is a boundary condition but it is actually an initial condition. You need boundary conditions at both ends of your domain. Your example code has alpha equal zero at both ends. Is this what you want?

[MT]

Dear Bill, I am sorry for the late reply, I had some trouble finding my own question back unfortunately. I had a look at the boundary conditions and I am not sure how to deal with that. Basically what I have is a column saturated with water (the alpha) which drains over time, so initially the whole column is saturated (alpha = 1) and in time it gets less since the fluid drains down and out of the column. In a discussion I had about this I was given the boundary condition alpha(x,t) = 0 for t>0 but I dont really understand this so I will have to ask again and then come back here, unless my explanation already told you something.

I appreciate your reply,

Martin