Solving PDE without BC with MATLAB

[smallphi]

I am trying to solve numerically the following PDE:

dF(x,t) / dt = some function of x and F(x,t) ONLY

where 0<x<5.

This equation does NOT need boundary conditions at x=0 and x=5 because each point in x evolves independently from the others (the equation doesn't contain spatial derivatives) and it is enough to know the initial condition F(x, t=0) to obtain the full solution F(x,t).

Unfortunately pdepe (the MATLAB numeric pde solver) asks for boundary conditions. I tried 0=0 but got an error. In principle I can evolve the points x=0 and x=5 and get the boundary conditions. I am trying to avoid that since for the actual problem I am trying to solve F is a column vector of a function, f(x,t) and its first and second spatial derivatives. Obtaining those at x=0 is tough since the function goes to zero so I will have to take a complicated limit which also depends on the behavior of f(x,t) around x=0 and f(x,t) is still unknown.

So is there a way to make pdepe forget about the boundary conditions?
Mathematica does it automatically.

 

[RedBranchKnight]

What is the precise form of the (linear??) PDE on the right-hand side of?

dF(x,t) / dt = some function of x and F(x,t) ONLY

If there are no spatial derivatives, then it is not a PDE, or am I missing something? Is it not just an ODE system? You also seem to contradict yourself later in your post (1st and 2nd order derivatives in x).
The variable x is just an extra parameter it seems?

 

[smallphi]

$$\frac{\partial R(x,t)}{\partial t} = \sqrt{\beta(x) + \alpha(x)/R}$$

$$P(x,t)=\frac{\partial R}{\partial x}$$

$$Q(x,t)=\frac{\partial P}{\partial x}$$P is the first derivative of R and Q is the second. Alpha and beta are some given functions. I could solve the first PDE as ODE in the time direction for each fixed x but that is inefficient, each call to ODE will take at least 0.01 sec and if I have 10 000 points that gives me about 1-2 minutes and that the optimistic expectation. This is not such a problem but then I will have to compute the derivatives P and Q from the R(x,t) data and I expect errors in that process, especially for the second derivative.

That's why I prefer to solve the whole problem for R and it's derivatives P and Q at once with the above system of PDEs hoping that the PDE solver will control the errors in P and Q. The initial condition is R(x,t=0)=R0(x). The initial conditions for P and Q are derived from that. Boundary conditions are not necessary for R(x,t) but I supplied anyways, I derived the BC for P and Q and was able to carry the computation in Mathematica. Matlab on the other hand is able to do the calculation for 50 points in x directions with obvious errors-oscillations around the actual solution. When I increased the points to 500 it started throwing errors that it can't meet the 'error tolerances close to t=0' it just stops the integration.

I wanted to do the problem in Matlab since its easier to program but will have to stick to Mathematica with its weird syntax but smarter routines.

 

[mercurial]

Hi. You're probably confusing yourself too much by wanting to think of this as a PDE. It really is just an ODE in disguise. That being said, Matlab should be able to solve it using ode45 or whatever solver provides a Runge-Kutta method. (Sorry, my Matlab is a bit rusty.) You just have to look at the problem right.

First question: how accurate do you need P and Q to be with with respect to that of R? Four decimal places? Ten? You can approximate both of those derivatives from your discrete values for R to order $$h^2$$ using a central differences. Let $$R = (R_1,\ldots,R_n)$$ be your approximation, at a given time step, on a grid $$x_1, \ldots, x_n$$ with mesh size $$h=x_{i+1}-x_i.$$ The respective approximations for $$\partial R/\partial x$$ and $$\partial^2 R/\partial x^2$$ are

$$P_i = \frac{R_{i+1} - R_{i-1}}{2h}$$

$$Q_i = \frac{R_{i+1} - 2R_i + R_{i-1}}{h^2},$$

both to accuracy $$O(h^2)$$. Thus if you need four decimal places of accuracy, chose a grid with mesh size $$h=10^{-2}$$, or about 500 points on your interval [0, 5].

For this your ODE would be just a first-order system with ~500 degrees of freedom with the nonlinear (time independent!) right hand side you've given. That's numerically tractable, unless you're integrating over a huge time scale.

If you need the same precision for P, Q at the boundary that you have in the interior, you will need to do something creative. Using the central difference approximations would require "phantom" points $$x_{-1} = -h$$ and $$x_{n+1} = 5+h$$ to make the formulae work at the boundary points.

I hope this gives you something to work with.

Jerry

 

[mercurial]

One more thing: you never said anything about the functions $$\alpha$$, $$\beta$$, or the initial condition $$R_0.$$ Presumably they're such that, say, the right hand side remains bounded away from zero and $$\sqrt{\beta + \alpha/R}$$ remains bounded for all time, preventing the solution from blowing up to $$\pm\infty.$$ How poorly behaved your RHS is may be the only thing that slows down numerical integration of the ODE system, especially if you're using an adaptive solver.

 

[smallphi]

Too much work with Matlab, I succeeded to get the function and the first derivative right but the second derivative shows strong oscillations. Also the calculation stops in the middle of the time interval cause 'the precision can't be satisfied'. On top of it its terribly slow ~ minutes.

Mathematica did it on the fly for 30 seconds, the function and the two derivatives included and they look smooth, no oscillations, without me bothering with automatic settings too much so I will have to stick to Mathematica.

It's a shocker to me that Mathematica(usually considered a symbolic program) has smarter and faster routines then Matlab(considered THE numeric program of choice).

 

[mercurial]

Hi Smallphi,

I'm glad you got Mathematica to produce an "answer," but do you have any way of checking it, aside from the fact that its first and second derivatives look smooth? I ask this because I've seen Mathematica produce wrong solutions to problems. Once I worked with a couple physicists who were using it to compute zero sets of systems of polynomial equations in the xy-plane. The more we examined the system, the harder it became to understand what Mathematica was producing, until we finally realized that it was wrong. In fact we found a mathematical proof that, with our particular system, Mathematica's alleged solutions were flat out impossible. Given that you seemed to have some confusion distinguishing an ODE from a PDE in your original post -- which may explain why one solver wouldn't handle it as a PDE -- you may want to verify your solution with another source. Remember: you're only looking for a solution of

$$\frac{\partial R}{\partial t} = \sqrt{\beta + \alpha/R},$$

which is an ODE because only the derivative w.r.t. one variable appears. The other variable x on which R depends is just a parameter here. The quantities P, Q depend only on the solution R, and not the other way around.

 

[smallphi]

1.I solved the system of R, P, Q equations as PDE system because that's what it is and because I wanted the PDE solver to control the numerical errors in P and Q.

2.I solved only the R equation as PDE cause it's faster than solving ODE at each x-point as I explained several times. And Mathematica is smart enough not to ask for 'boundary conditions' for that PDE while Matlab refuses to solve unless I specify some fictitious BC.

3.Moreover I have analytic expressions for R, P, Q. I don't use them afterwards because later I have to solve other PDE's involving R, P, Q and it's faster to calculate R, P, Q once and for all as interpolation functions than calling the analytic routines every time. Plus the analytic routines have problems at certain points, have to take limit when x->0 etc.

The results of all 3 methods coincide within numerical errors so I have no doubts Mathematica is right. Matlab concurs with Mathematica for R and its first x-derivative P, after one day of playing with settings and still can't integrate it for the full time interval. The second derivative in Matlab was a mess and I don't have a week to play with settings to get it right. Why bother when Mathematica solves the whole problem as one system of PDE's for 30 secs.

I've noticed many times the Matlab routines are not as smart as Mathematica's handling various situations automatically. Like for example, in Matlab there's no such thing as 'integration to infinity' which Mathematica handles automatically. That is a sad conclusion since programming in Matlab is way more user-friendly than Mathematica. I really think that the people programming the Matlab routines should write smarter, more automatic routines than leaving every annoying detail to the user. And please don't try to make it look more 'object-oriented' like Mathematica, that complicates the programming and instead of focusing on the problem at hand, the user has to bother with programming syntax like in Mathematica. One hour reading of help to figure out how a specific command should be written ... Nobody needs that and the newer versions of Matlab get worse in that aspect.

When I asked the question on the Matlab forum some fool thought I was spamming for Mathematica LOL I use all three programs Mathematica, Maple and Matlab. For symbolic computation I prefer Maple, for numeric Matlab in general unless its problematic like in this case and I have to go with Mathematica. I wish they wrote Mathematica in a simple language like Maple instead of those stupid Lisp-like constructs.