How to formulate boundary conditions for a PDE system?

I have the following PDE system steaming from Flash Photolysis:

I formulated the PDEs in the form reqired by pdepe and coded pdefun as follows:

   function [C,F,S] = photopde (x,t,U,DUDX)
        C = [0;1];
        F = [0;0];
        S = [epsilon*U(1)*U(2);phi*epsilon*U(1)*U(2)];
   end

The initial coditions are:

which are encoded into an icfun:

function U = photoic (x)
      U = [I0*exp(-C0*epsilon*phi*x);C0];
end

However I have only 2 boundary conditions at the left side at t = 0 and none on the right side:

How do I code corresponding bcfun? How do I define PR and QR outputs from bcfun? Thanks.

Answer 1

[Torsten]

Don't use "pdepe" to solve this system.

Discretize dI/dx in space. This will result in algebraic equations for I in each spatial grid point.

Then use ODE15S to solve the second ordinary differential equation for C in each of these spatial grid points.

Best wishes

Torsten.

[Mikael Hakman]

Torsten,

I don't know how discretize dI/dx because I don't know C(x,t) - this and I(x,t) is what MATLAB should give me, numerically. For example at x=0.5 I get:

dI/dx(0.5,t) = ε*I(0.5,t)*C(0.5,t)

What I'm supposed to do with this? If you know how to discretize then please show me. Note that I have the actual analytical solution but that is beyond the point. Thanks.

[Torsten]

You can use this function for the ODE integrator ODE15S to calculate C:

function dCdt = myfun(t,y)

C=y;
I(1)=I0;
for k=2,length(C)
  I(k) = I(k-1)+deltax*(-epsilon*I(k-1)*C(k-1));
end 
dCdt(1,1)=C0*(-J0*phi*epsilon)*exp(-J0*phi*epsilon*t);
for k=2,length(C)
  dCdt(k,1)=-phi*epsilon*I(k)*C(k);
end

Best wishes

Torsten.

[Mikael Hakman]

What you are suggesting is that I should integrate I(x,t) by "hand" using the simplest of all integration methods - the Euler method. Then why not do the same for C(x,t) if that is what was required. But then we can do that in C++ or Java as easy as in MATLAB, don't we?

[Torsten]

It's the method-of-lines approach - the same approach "pdepe" would use if it was able to solve your equations.

My advise is to start simple. You may try a 2nd order scheme in space for I later, but then it must be an upwind scheme for stability reasons. Centered differences will lead to oscillations, I guess.

Best wishes

Torsten.

Answer 2

[Bill Greene]

pdepe is specifically designed for PDEs that are second-order in the spatial direction. That is why it requires you to specify boundary conditions at both ends.

But that doesn't mean it won't solve first-order PDEs and, since you already have most of the code implemented, I suggest you try it. Since your flux is zero everywhere, a boundary condition you can use at the right end to satisfy pdepe is qr=1, pr=0.

Sometimes when solving first-order PDEs with pdepe the solution shows some spurious oscillations. One simple way to deal with this problem is to add a small amount of "artificial diffusion", e.g.

F = 1e-6*[1;1];

where the factor can be adjusted so that it damps out the oscillations without having a large effect on the solution.

[Mikael Hakman]

I tried this method before I asked the question here. I get the following picture from MATLAB:

The analytical solution gives:

I'm attaching full MATLAB source here. Thanks.

[Hitesh Saini]

Hi Bill

How can we apply boundary conditions while solving PDE using line of approach method?

Is it possible?