pde1dm compared to pdepe

pde1dm is an extended version of pdepe and it can solve equations pdepe can't solve. I wonder about the differences between these two algorithms in the spatial discretization?

[feynman feynman]

It's said that pdppe and pde1dm are both based on finite elements. But sometimes when the spatial discretization is too coarse, both will have difficulty converging and thus report errors. This doesn't seem to me to me a feature of finite elements because even if the mesh is very coarse in finite elements there is no nothing to stop the program from running and cause reporting errors. Could anyone help thank you!

[Torsten]

If the mesh is too coarse, solutions can diverge and thus make the ODE solver to quit.

I don't know where you found that "even if the mesh is very coarse in finite elements there is no nothing to stop the program from running and cause reporting errors."

[feynman feynman]

With my limited experience with finite elements I think there is no such requirement on the spatial step length as in that in finite difference time integration for its stability?

[Torsten]

How should the finite element method produce a result converging to a senseful solution if you use e.g. (exaggerated, I admit), three mesh points for a distance of 100 km ? I can imagine that such a problem could produce Inf or an oscillating solution in the center point during time integration.

And there is no result I remember that says "the finite element method is stable per se". This can't be true because there is no such thing as "the" finite element method.

[feynman feynman]

right, thank you!

Answer 1

[Bill Greene]

There are several differences in the pde1dm spatial discretization compared with pdepe.

For the m=0 case (Cartesian), both pde1dm and pdepe use the standard, linear shape functions. But, by default, pde1dm evaluates the residual using a two-point Gaussian integration rule while pdepe evaluates only at the center of the element. Often this doesn't significantly affect the solution but, for some problems, this allows pde1dm to obtain a converged solution when pdepe cannot.

Of course, by default, this makes pde1dm somewhat slower than pdepe but use of the "Vectorized" option with pde1dm dramatically reduces the computational time.

In evaluating the so-called mass terms (those arising from the c-coefficients), pdepe essentially uses what is often referred to as a "lumped" approach that results in a diagonal mass matrix. pde1dm evaluates these terms using exactly the same approach as the f and s terms ("consistent" approach).

For cylindrical (m=1) and spherical (m=2) geometries, pdepe uses special shape functions, as described in the Skeel and Berzins paper referenced in the pdepe documentation. According to the paper, this improves the convergence near the singular point (x=0). pde1dm does not use these special shape functions.

The above describes the current implementation of pde1dm but may change in the future. Also, pde1dm has some undocumented options that cause pde1dm to produce a discretization essentially the same as pdepe for the m=0 case. But these options may change or be eliminated in future versions; hence they are undocumented for now.

[feynman feynman]

Thanks very much for the detailed answer. If both these algorithms are based on finite elements, why are sometimes initial conditions automatically adjusted through iterations so that u(t=0) isn't the same as what the user assigns?

When I was trying to use the vectorized option, it seems pde1dm(m,@pde,@ic,@bc,x,t,options) returns errors, while pdepe(m,@pde,@ic,@bc,x,t,options) doesn't. That is, options isn't recognized or authorized by pde1dm in the above form, and perhaps it needs previous steps before this line?

[Bill Greene]

why are sometimes initial conditions automatically adjusted through iterations so that u(t=0) isn't the same as what the user assigns?

If you include an elliptic equation (c=0) in your PDE system the discretization produces a system of DAE rather than ODE. In general, to solve a DAE using ode15s or ode15i, you must have a "consistent" set of initial conditions at t=0; that is they must satisfy the DAE system. So the initial conditions are adjusted automatically to insure this. If you want to know more, look at the Matlab documentation for these ODE solvers and the referenced papers by Shampine.

When I was trying to use the vectorized option, it seems pde1dm(m,@pde,@ic,@bc,x,t,options) returns errors, while pdepe(m,@pde,@ic,@bc,x,t,options) doesn't. That is, options isn't recognized or authorized by pde1dm in the above form, and perhaps it needs previous steps before this line?

Are you using odeset to set the options argument? pdepe allows the options argument to be set by the odeset function. odeset sets many options for all the Matlab ODE solvers but pdepe handles only a few of these and ignores the rest. pde1dm allows only the options that it actually supports. For example if you want to change RelTol, just define options as

options.RelTol=1.e-4;

[feynman feynman]

thank you so much! I basically just:

options=odeset('RelTol',1e-3)

sol=pde1dm(0,@pde,@ic,@bc,x,t,options)

but it shows an error: "'Events' is not a recognized parameter. For a list of valid name-value pair arguments, see the documentation for this function." This won't happen when pde1dm is replaced with pdepe.

[Torsten]

I interpreted the last response that you have to use

options.RelTol=1.e-4;
sol=pde1dm(0,@pde,@ic,@bc,x,t,options)

instead of

options=odeset('RelTol',1e-3)
sol=pde1dm(0,@pde,@ic,@bc,x,t,options)

[feynman feynman]

it worked, thanks very much! I don't get it though, because for pde1dm, 'options' wasn't even defined before 'options.RelTol=1.e-4;' I guess 'odeset' doesn't play with pde1dm.

[Torsten]

For "pdepe" and "pde1dm", "options" is just a structure variable. For "pdepe", it can be created via "odeset", for "pde1dm", you should create it on your own. If you'd create it via "odeset" for "pde1dm", some of the structure components below were most probably missing in "pde1dm" (e.g. "Events") and could not be interpreted.

odeset

          AbsTol: [ positive scalar or vector {1e-6} ]
          RelTol: [ positive scalar {1e-3} ]
     NormControl: [ on | {off} ]
     NonNegative: [ vector of integers ]
       OutputFcn: [ function_handle ]
       OutputSel: [ vector of integers ]
          Refine: [ positive integer ]
           Stats: [ on | {off} ]
     InitialStep: [ positive scalar ]
         MaxStep: [ positive scalar ]
             BDF: [ on | {off} ]
        MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]
        Jacobian: [ matrix | function_handle ]
        JPattern: [ sparse matrix ]
      Vectorized: [ on | {off} ]
            Mass: [ matrix | function_handle ]
MStateDependence: [ none | {weak} | strong ]
       MvPattern: [ sparse matrix ]
    MassSingular: [ yes | no | {maybe} ]
    InitialSlope: [ vector ]
          Events: [ function_handle ]

[feynman feynman]

oh so well explained thanks!

[feynman feynman]

I asked these in separate ones but was suggested to ask in this same request

pde1dm with coupled ODEs - MATLAB Answers - MATLAB Central (mathworks.com)

pdepe/pde1dm's vectorized option - MATLAB Answers - MATLAB Central (mathworks.com)

[feynman feynman]

@Bill Greene

Dear Dr Greene, normally pde1dm works properly but very often it reports the following errors when closing matlab and reopening it. I wonder if there's anything wrong?

Not enough input arguments.

Error in
PDE1dImpl>@(t,y,yp)self.calcResidual(t,y,yp,depVarClassType)
(line 158)
icf = @(t,y,yp)
self.calcResidual(t,y,yp,depVarClassType);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode15i (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0,
options, varargin);

Error in PDE1dImpl/solveTransient (line 253)
[outTimes,u]=ode15i(icf, self.tspan, y0, opts);

Error in pde1dm (line 125)
u = pdeImpl.solveTransient(odeOpts);

Error in mypde (line 47)
sol=pde1dm(m,@pdefun,@icfun,@bcfun,x,t,options);

[Bill Greene]

Post your code that causes this error.

I have no idea what you mean by "when closing matlab and reopening it" so, if sopmething special needs to be done to causes this error, I'll need more details.

[feynman feynman]

There's nothing special about my code, just your pde1dm package folder with the following main code. When this error occurs, simply closing and reopening matlab will cure it. Then closing and reopening matlab again might cause this error to reappear.

function heat
x=linspace(-pi,pi);t=linspace(0,10);
m=0;sol=pde1dm(m,@pde,@ic,@bc,x,t)
u=sol(:,:,1);
surf(x,t,u)
    function [c,f,s]=pde(x,t,u,ux)
        c=1;f=ux;s=0;
    end
    function U0=ic(x)
        U0=cos(x);
    end
    function [pl,ql,pr,qr]=bc(xl,ul,xr,ur,t)
        pl=ul+1;pr=0;ql=0;qr=1;
    end
end

[Bill Greene]

I did the following:

1. Downloaded a fresh copy of pde1dm from github.
2. Created a new file, heat.m, in the pde1dm directory.
3. Pasted your code into this file.
4. Ran your code in heat.m several times with no error.
5. Restarted matlab.
6. Ran your code several more times with no error

I've looked at the lines of code pointed to by your error messages and I simply have no clue what could be causing the problem.

I am running matlab R2019b, by the way.

[feynman feynman]

Thank you for your time! I use r2021a, is it unstable?

[feynman feynman]

I'm calling different functions in different folders via

switch solver
    case 'solver1'
        addpath 'D:\matlab\project1\';
        solver1
    case 'solver2'
        addpath 'D:\matlab\project2\';
        solver2

Could this frequently using addpath cause any problems?

[feynman feynman]

@Bill GreeneIn addition, is ode15i in your pde1dm irreplacable with any other ode solvers like the following?

[outTimes,u]=ode45(icf, self.tspan, y0, opts);

% [outTimes,u]=ode15i(icf, self.tspan, y0, yp0, opts);

[Torsten]

ode15i is an implicit ODE integrator and has a substantially different interface compared to the other solvers. It cannot be replaced by another integrator from the ode suite.

[feynman feynman]

I saw there was an if statement about whether to use ode15i and the code requires that this statement be always true. I wonder if one could add functionality so that other ode solvers can be put into the else statement.

Answer 2

[Torsten]

I wonder about the differences between these two algorithms in the spatial discretization?

No difference. pde1dm has the option to add ordinary differential equations to the system of partial differential equations and in the boundary condition part, but this doesn't change the numerical kernel.

[feynman feynman]

Thank you. However when the same equation is solved the two solvers give quantitative slightly different results, that's why I doubted. If no difference, how does pde1dm manage to be more flexible in incorporating ODEs?

[Torsten]

You should contact the author of the code - Bill Greene.