PDEPE for simple coupled ODE-PDE system?

Reposting (just once) to see if anyone has any insight. I solved the problem using the "method of lines" technique but was really curious about my original question and the interesting behavior I observed.

Hi all, I have a system of coupled reaction-diffusion ODEs and PDEs that I wish to use pdepe to solve, and have been running into unexpected outputs. To check, I used a simpler "test" system which could be solved in an exact fashion, purely to test the output of PDEPE.

It is simply:

 function [cc,ff,ss]=brpde(x,t,u,dstatedx)
 mu=0.02;
 u1=u(1);
 u2=u(2);
     if (x >5)
        du1 = -0.07*u1;
        else
            du1=0;
        end
        if (x > 5)
        du2=-0.07*u2;
        else
            du2=0;
        end
 ss=[du1; 
    du2];
 ff=[0;
    mu*dstatedx(2);];
 cc=[1;
    1;];
 return
I have set the boundary conditions to impose no flux and the left and right hand boundaries (this should work, I believe, even for the ODE variable u1 where flux is zero throughout the slab).

function [pa,qa,pb,qb]=brbc(xa,ua,xb,ub,t) vecone=ones(2,1); veczero=zeros(2,1); pa =veczero; qa = vecone; pb = veczero; qb = vecone; return And suitable initial conditions:

function u0=bric(x) u0=[10; 10;]; return The output for u2 looks correct as I would expect it. However, although the exact solution is that u1 (the variable with no flux) should asymptotically approach zero for x>5 and stay at 10 for x<5, at the 'interface', ie near x=5, the solutions transiently dip below zero (even when the mesh is made finer or the error tolerances are made more stringent). It's confusing behavior, and I am wondering whether it's since pdepe is not suited for mixed ODE-PDE systems. When I code this simple system into a discretized forward Euler, of course the output is as I would predict it.

Thought I would reach out to the experience of the community. Anyone have suggestions? Or is there an obvious error in the above? (My real application has coupling between the ODEs and PDEs, but if it doesn't work for this test case, I am leaning toward hard-coding it).

Thanks in advance! VI

Answer 1

[Torsten]

I guess the behaviour stems from the discontinuity of u1 at x=5.

Try

du1 = -0.07*u1

instead of

if (x >5)

du1 = -0.07*u1;

else

du1=0;

end

in your code to see if my guess is correct.

Best wishes

Torsten.

[Vivek]

This sure worked, thanks for the suggestion on troubleshooting!

For problems in which there is spatial heterogeneity in the reaction terms (or say, heat generation in a wire of varying radial properties), or such as the discontinuity I introduced, is there a common workaround you can think of that would permit use of pdepe? (I am thinking of something like the 'Events' flag in odeset, except defining critical radii or dimensions at which discontinuities may occur).

The method of lines technique I tried seemed to produce expected solutions even with simple forward Euler for spatial integration in this "test problem", so perhaps a hand-crafted solution is better. Not sure why that would handle discontinuities better, but perhaps the answer is in how the two techniques differ in handling spatial steepness?

Regardless, thanks for the answer and happy new year. VI