Compare PDESolve

JohnRudnicki · ‎Oct 10, 2009

I have posted about this before, but I think this is the end. The attached worksheet compares solutions for a system of two, linear diffusion-like pdes by pdesolve, an explicit finite difference method, an implicit Crank-Nicolson method and the analytic (exact) solution. The exact solution decays exponentially (for the wave number used), but all three numerical methods diverge (grow large) after some time. For the same number of spatial points, and roughly the same time increments, the explicit FD does better than pdesolve. The Crank-Nicolson method does the best with much larger time steps, but also diverges.

The divergence of the numerical solution is not surprising (in retrospect) since the actual solution is so small by this time that it is likely on the order of the numerical errors. But the unwary might conclude that it is the pde's that are unstable and they are not.

Any suggestions on improved programming are welcome!!

AlvaroDíaz · ‎Oct 10, 2009

On 10/10/2009 6:10:13 PM, jwrudn wrote:
>I have posted about this ...

I see an extra point in the link to the mcd file, maybe this is the reason for that I can't download it.

Alvaro.

ptc-1368288 · ‎Oct 11, 2009

On 10/10/2009 7:09:59 PM, adiaz wrote:
>On 10/10/2009 6:10:13 PM, jwrudn wrote:
>>I have posted about this ...
>
>I see an extra point in the link to the
>mcd file, maybe this is the reason for
>that I can't download it.
>
>Alvaro.
__________________________

The collab must remove the extra (.) in the sheet name.

jmG

JohnRudnicki · ‎Oct 11, 2009

Yes, sorry. It was the end of a long day. Re-posting without the extra dot.

ptc-1368288 · ‎Oct 12, 2009

On 10/11/2009 2:52:25 PM, jwrudn wrote:
>Yes, sorry. It was the end of
>a long day. Re-posting without
>the extra dot.
_______________________________

You don't arrive at the PDE's from analytical development and how they should appear in the solver. The work sheet must be reconstructed for versions lower than 14 (?), because of the inclusion of so many definitions in the solver. Why do you have to ceil/floor in the first place ? just out of the blue. In the past (the other work sheets), I mentioned the "spatial had to be 1, maybe 0.5... and because you have two DE's the solver will not reply "argument outside of range". Past the tutorial (a bit dense for the first time reader), then there are many pages missing ( graphics, infinitesimal analysis, uniqueness/existence ...).

Just comments suggesting that you have to compose the work sheet from blank. What you have there is neither a piece of maths, nor a piece of Engineering.

jmG

ptc-1368288 · ‎Oct 12, 2009

You should try understanding what you are trying to solve for. Make the attached work from other inputs in the system and don't plot, export as soon as posible, plots are for the birds if you can't export from them. What has Fourier to do in there ?

jmG

JohnRudnicki · ‎Oct 13, 2009

On 10/12/2009 2:58:31 AM, jmG wrote:
>You should try understanding
>what you are trying to solve
>for.

I do know what I am solving for: pressure and temperature, denoted by uppercase psi and phi in pdesolve. The �analytic� solution is correct. I don�t derive it here, but you can check it, if you want to. Look for solutions of the form exp(pt)cos(2 pi n Y). Solution of the quadratic equation for p is given in the worksheet. n must be integer to satisfy the bc. Your plot of pplus and pminus for n from -10�10 shows that the real part of plus and pminus is positive as it should be for n less that the critical value, 10.9 for the values used.

The finite difference approximations also seem to be correct (though probably can be programmed more efficiently) since they track the analytical and pdesolve solutions well at least if time is not to large and n is not too small.

(worksheet is not set to calculate automatically since, depending on the number of points, some of the calculations run for a bit).

>Make the attached work
>from other inputs in the
>system and don't plot, export
>as soon as posible, plots are
>for the birds if you can't

Not sure what you mean here.(I manage to use MC but I am certainly not an expert) Cant you just save the worksheet and alter the inputs.

>export from them. What has
>Fourier to do in there ?

By �Fourier� I simply meant a cos or sin with a single wavelength, maybe not the best term.

jwr

JohnRudnicki · ‎Oct 13, 2009

On 10/12/2009 12:44:53 AM, jmG wrote:
>On 10/11/2009 2:52:25 PM, jwrudn wrote:
>>Yes, sorry. It was the end of
>>a long day. Re-posting without
>>the extra dot.
>_______________________________
>
>You don't arrive at the PDE's from
>analytical development and how they
>should appear in the solver.

The worksheet was not meant as a tutorial (not that I could give one on numerical solutions of pdes). I have not gone into the physical derivation of these equations because I did not think it was relevant here.

>The work
>sheet must be reconstructed for versions
>lower than 14 (?), because of the
>inclusion of so many definitions in the
>solver.

I haven�t used MC11 in years and would not know how to rewrite the worksheet for this version.

>Why do you have to ceil/floor in
>the first place ? just out of the blue.

n must be integer to satisfy the bc: floor and ceil pick out integers just above and below the critical n (which is generally not integer).

I dont know what you mean by..

>In the past (the other work sheets), I
>mentioned the "spatial had to be 1,
>maybe 0.5... and because you have two
>DE's the solver will not reply "argument
>outside of range".

Spatial range is set by the physical problem and nondimensionalized to 0�1

> Past the tutorial (a
>bit dense for the first time reader),
>then there are many pages missing (
>graphics, infinitesimal analysis,
>uniqueness/existence ...).
>
>Just comments suggesting that you have
>to compose the work sheet from blank.
>What you have there is neither a piece
>of maths, nor a piece of Engineering.
>
>jmG

jwr

ptc-1368288 · ‎Oct 13, 2009

On 10/13/2009 1:22:53 PM, jwrudn wrote:
...
>The worksheet was not meant as a
>tutorial (not that I could give one on
>numerical solutions of pdes). I have not
>gone into the physical derivation of
>these equations because I did not think
>it was relevant here.
>
...
>jwr
___________________
"The worksheet was not meant as a tutorial"
...then : the work sheet was not meant so.

If you don't derive the "physical infinitesimal elemental differences", how can you decide about the differential representation ? then the composition of the system ? then its representation, then the solution by either method(s) ??? Start by drawing a picture of the system and the difference elements, then the remaining of the analysis. You can't tell in advance of the "existence, bounds, uniqueness", whether the system has or nor a numerical solution. The way I recall floor/ceiling, it governs the "spatial=1", it's that spatial that you have to mesh and scale to the physical dimension, not the way you are saying. Your system is only valid within floor/ceiling, you can view as trig in rad. But your system is specific to user. You can think also floor/ceiling like Chebyshev range.

jmG

TomGutman · ‎Oct 11, 2009

I will disagree with your charactization of the PDE. I consider them unstable, generally going to infinity.

But the stability is dependent on the initial conditions. For pure harmonic initial conditions the equations are stable for sufficiently high frequency (actually decays to zero).

But for any actual use, whether real or numeric, the instability will dominate. Reality is never so nice as to be strictly bandwidth limited, and I would expect some low frequency components in any real system. With the unstable PDEs there is no such thing as an ignoribly small low frequency component. Numerically you have a similar situation. Even if you attempt a pure high frequency initial condition, numerical errors, both in the original representation and introduced at each step will inevitably introduce low frequency components and these will eventually dominate the solution.
__________________
� � � � Tom Gutman

JohnRudnicki · ‎Oct 11, 2009

Tom, thanks for your comment. I agree, you are correct. When I said that the equations are stable, I meant that the exact solution of the equations is stable (decays exponentially) for the given initial conditions used in the worksheet. As you say, stability does depend on the initial conditions. For a single Fourier mode, they are only stable for sufficiently short wavelengths (low frequencies). Your discussion clarified what I intuitively thought might be going on.

More generally, however, my concern is when there is no exact (analytic solution) how do you tell whether instability is numeric or due to the equations themselves (without knowing a lot more about numerical analysis,error propagation, etc. than I do)? Seems like you could have the first w/o the second. Probably not a question that is easily answered but for me, the question is not entirely academic. I am interested in a set of nonlinear pdes of which the version in the worksheet is linearized.

AlvaroDíaz · ‎Oct 12, 2009

On 10/11/2009 7:50:32 PM, jwrudn wrote:
>Tom, thanks for your comment. I agree, you are correct. When I said that the equations are stable, I meant that the exact solution of the equations is stable (decays exponentially) for the given initial conditions used in the worksheet.

This isn't the concept of stability. If a magnitude grows without limit (eventually exponentially) have a physical sense: something go to broken, and this not have relation with the stability of the solutions. In your case: if the pressure grows too much, then sure the material under this pressure go to have a fracture. And then your model must to change, because now is not true or good for the new situation. Same thing with the temperature, let say that the material could burn or change the state, or something that make the model inadequate and need to set to another model.

Also this not meaning that the model is not useful: only says that you need to check the material properties to know the time limit that could be exposed to this situation before it go to be broken.

For a rigourous definition for the 'stability in the future as Lyapunov' you can see

http://www.fing.edu.uy/~eleonora/dvi/cual4.pdf

in spanish, but for odes, not pdes.

Stability is usually the study of small variations of the parameters in the diff equations, and is usual taking initials conditions as parameters, but eventually could be anything. For example, c_thermical or c_hydraulic could be the parameters to study the stability. In your case seems that the stability that you want to study is precisely this 'stability in the future', this is, for t big (4 or 5 seconds are enough big with the others parameters fixed). So, your parameter is the time, and the 'small variations' about a fixed point is the grow of the time without limit. As you can see in the paper, particular solutions (which depends on initials conditions) must to be taken in a 'maximal interval', to ensure that the solution is stable also for small variations around the initials conditions. So, this stability must to study first.

In your case, for N>Ncrit=10.972 what you have is not stability, only solutions that not diverges to inf as t diverges to inf. This meaning only that your material not burns or fissure, but not implies that the analytical solutions are stables. Notice that for N=6, for example, analytic and pdsolve solutions are the same.

For numerical solutions the situation is more complicated: you need that the analytic solution must to be stable, but also needs to add numerical criterion to the 'stability' of the numerical solutions, but this depends on the particular algorithm, and mathcad don't give information about their numerical algorithms to do this.

Regards. Alvaro.

TomGutman · ‎Oct 12, 2009

It's even a bit more complicated. Even when the analytic solution is stable, a numeric solution mught not be. But there are two distinct causes for this. First is the numeric algorithm. Typically the errors implied by the algorithm are the dominant errors and the behaviour of the algorithm is the important issue. But there is also the numeric calculation that also introduces errors. Common cause is where the algorithm requires the difference between two quantities that are very close together. Mathematically, that is not a problem. Numerically, the limited precision of the number representation can lose all precision. So you have to worry first as to whether the numeric algorithm works mathematically, and then whether it has problems with the limited computer arithmetic that is actually used for the calculations.
__________________
� � � � Tom Gutman

ptc-1368288 · ‎Oct 12, 2009

>It's even a bit more complicated<<br> __________________________________

Surely: if the two "elemental solutions" are as given, then where is the XFR function as the resulting differential PDE system ? Answer: nowhere in the work sheet.

jmG

JohnRudnicki · ‎Oct 13, 2009

Alvaro: Thanks for your comment. I agree with you though your notion of �stability� seems more sophisticated than mine. When solutions grow exponentially (as they do for initial conditions with large enough wavelengths, small enough N�s), there is no question of �numerical stability�. This growth does have physical significance, though I have not gone into it here and it is a result of linearizing the original system. In the actual nonlinear system this growth is, I think, and, as you imply, mitigated by other physical factors in the nonlinear system (though I still need to solve this which is why I wanted to study the linear system to get some experience with numerical solution.)

When the initial conditions have wavelengths that are short enough (large enough N�s) the exact solution decays exponentially, but near the critical N your comments about stability are pertinent and I think related to Tom�s point, if I understand it. Numerical approximations will cause the (numerical) solution to become unstable. And, yes, the system is unstable for initial conditions that contain all wavelengths

jwr

ptc-1368288 · ‎Oct 14, 2009

On 10/10/2009 6:10:13 PM, jwrudn wrote:
>I have posted about this before,
==> YES, under another sheet name... now untraceable !

>but I think this is the end.
OK: the end.
Just remember that some systems might "survive non sense". In the attached Cauchy predator/prey, the system survives with only 1 (one) wolf. I can take it, but not very honorable for the "Wolves" as they would have to all be hermaphrodite. For sure, the 1(one) wolf system is not Canadian !

jmG