Proof of Harten’s Lemma re the Convergence of TVNI Finite-Difference Schemes

Don Warrington Fluid Mechanics, Mathematics Leave a comment

An academic paper with a rather unusual history is that of the Israeli mathematician Ami Harten’s “High Resolution Schemes for Hyperbolic Conservation Laws.”  First published in the Journal of Computation Physics in 1983, it was republished in 1997 in the same journal, and is often cited with the later date.

At the time of republishing Peter Lax, who earlier had saved the computer from the hippie radicals, made the following statement about this paper in an introduction:

This paper was a landmark; it introduced a new design principle—total variation diminishing schemes—that led, in Harten’s hands, and subsequently in the hands of others, to an efficient, robust, highly accurate class of schemes for shock capturing free of oscillations. The citation index lists 429 references to it, not only in journals of numerical analysis and computational fluid dynamics, but also in journals devoted to mechanical engineering, astronautics, astrophysics, geophysics, nuclear science and technology, spacecraft and rockets, plasma physics, sound and vibration, aerothermodynamics, hydraulics, turbo and jet engines, and computer vision and imaging.

One point in the paper was a lemma concerning the convergence of TVNI (total variation nonincreasing) finite difference schemes.  Concerning the name of these schemes, Lax points out the following:

Harten originally called his schemes variation diminishing, abbreviated TVD; when Osher pointed out the usual meaning of these initials, the name was switched to total variation nonincreasing (TVNI), but was eventually settled on the more euphonious TVD.

The following is an expansion of Harten’s proof of the lemma.

Schemes which are total variation nonincreasing (TVNI) can be characterized as follows:
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|\leq\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|

We can thus define
TV\left(u^{n}\right)=\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|

TV\left(u^{n+1}\right)=\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|

and substituting
TV\left(u^{n+1}\right)\leq TV\left(u^{n}\right)

Consider the general expression
u_{{j,n+1}}=u_{{j,n}}-C_{{j-1,n}}\left(u_{{j,n}}-u_{{j-1,n}}\right)+D_{{j,n}}\left(u_{{j+1,n}}-u_{{j,n}}\right)

where
C_{{j-1,n}}\geq 0

D_{{j,n}}\geq 0 C_{{j-1,n}}+D_{{j,n}}\leq 1

We should observe that our ultimate goal is to sum these values from negative infinity to positive infinity; thus, we can shift the index at will. The inequalities will still hold but the specific location in space may change. It is also worth noting that the coefficients may themselves change at different points in space.

Let us consider the next spatial step, to wit
u_{{j+1,n+1}}=u_{{j+1,n}}-C_{{j,n}}\left(u_{{j+1,n}}-u_{{j,n}}\right)+D_{{j+1,n}}\left(u_{{j+2,n}}-u_{{j+1,n}}\right)

Subtracting the previous spatial step from this yields
u_{{j+1,n+1}}-u_{{j,n+1}}=u_{{j+1,n}}-C_{{j,n}}\left(u_{{j+1,n}}-u_{{j,n}}\right)+D_{{j+1,n}}\left(u_{{j+2,n}}-u_{{j+1,n}}\right)-u_{{j,n}}+C_{{j-1,n}}\left(u_{{j,n}}-u_{{j-1,n}}\right)-D_{{j,n}}\left(u_{{j+1,n}}-u_{{j,n}}\right)

Some rearranging yields
u_{{j+1,n+1}}-u_{{j,n+1}}=\left(u_{{j+1,n}}-u_{{j,n}}\right)\left(1-D_{{j,n}}-C_{{j,n}}\right)+C_{{j-1,n}}\left(u_{{j,n}}-u_{{j-1,n}}\right)+D_{{j+1,n}}\left(u_{{j+2,n}}-u_{{j+1,n}}\right)

Taking the absolute value of both sides, we have
\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|=\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\left(1-D_{{j,n}}-C_{{j,n}}\right)+C_{{j-1,n}}\left(u_{{j,n}}-u_{{j-1,n}}\right)+D_{{j+1,n}}\left(u_{{j+2,n}}-u_{{j+1,n}}\right)\right|

At this point we observe that
C_{{j-1,n}}\geq 0
D_{{j+1,n}}\geq 0
D_{{j,n}}+C_{{j,n}}\leq 1
1-D_{{j,n}}-C_{{j,n}}\geq 0

We can thus limit the absolute values and write the expression as follows:
\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|\leq\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|\left(1-D_{{j,n}}-C_{{j,n}}\right)+\left|\left(u_{{j,n}}-u_{{j-1,n}}\right)\right|C_{{j-1,n}}+\left|\left(u_{{j+2,n}}-u_{{j+1,n}}\right)\right|D_{{j+1,n}}

Taking the summation for both sides,
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|\leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|\left(1-D_{{j,n}}-C_{{j,n}}\right)+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j,n}}-u_{{j-1,n}}\right)\right|C_{{j-1,n}}+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+2,n}}-u_{{j+1,n}}\right)\right|D_{{j+1,n}}

Since, as we observed before, we can shift the indices (as the “centre” of the system is arbitrary with infinite boundaries) we can rewrite the above as follows
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right| \leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|\left(1-D_{{j,n}}-C_{{j,n}}\right)+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|C_{{j,n}}+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|D_{{j,n}}

in which case
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right| \leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|

Substituting, we have at last
TV\left(u^{n+1}\right)\leq TV\left(u^{n}\right)

Citation: Harten, A. (1997) “High Resolution Schemes for Hyperbolic Conservation Laws.” Journal of Computational Physics, Vol. 135, pp. 260-278.

Advertisements

Leonhard Euler on the Creator and Mathematics

Don Warrington Mathematics Leave a comment

Quoted in S. Timoshenko’s History of Strength of Materials:

Since the fabric of the universe is most perfect, and is the work of a most wise Creator, nothing whatsoever takes place in the universe in which some relation of maximum and minimum does not appear.  Wherefore there is absolutely no doubt that every effect in the universe can be explained as satisfactorily from final causes, by the aid of the method of maxima and minima, as it can from the effective causes themselves…Therefore two methods of studying effects in nature lie open to us, one by means of effective causes, which is commonly called the direct method, the other by means of final causes…One ought to make a special effort to see that both ways of approach to the solution of the problem be laid open; for this not only is one solution greatly strengthened by the other, but, more than that, from the agreement between the two solutions we secure the very highest satisfaction.

Very few people on the earth have contributed to the basics of mathematics and mechanics than Euler.  Living in the “Age of Reason” Euler pushed science forward while remaining a Christian all of his life.

Direct Derivation of the Equation of Motion for an Undamped Oscillating System in Phase Angle Form

Don Warrington Mathematics, Solid Mechanics 1 Comment

The equations of motion for linear vibrating systems are well known and widely used in both mechanical and electrical devices. However, when students are introduced to these, they are frequently presented with solutions which are either essentially underived or inadequately so.

This brief presentation will attempt to address this deficiency and hopefully show the derivation of the equation of motion for an undamped oscillating system in a more rigourous way.

Consider a simple spring/mass system without a forcing function. The equation of motion can be expressed as

1

where x(t) is displacement as a function of time, m is the mass of the system, and k is the spring constant. The negative sign on the right hand side of the equation is not an accident, as the spring force always opposes the motion of the mass, and is the result of using a mechanical engineers’ “free body diagram” method to develop the equation.

Solutions to this equation generally run in two forms. The first is a sum of sines and cosines:

2

But it’s more common to see it in the form of

3

The latter is simpler and easier to apply; however, it is seldom derived as much as assumed. So how can it be obtained from the original equation?

Let us begin by considering the original differential equation. With its constant coefficients, the most straightforward solution would be a solution where the derivative (and we, of course, would derive it twice) would be itself. This is the case where the function is exponential, so let us assume the equation to be in the form of

4

(I had an interesting fluid mechanics/heat transfer teacher who would say about this step that “you just write the answer down,” which we as his students found exasperating, but this method minimises that.)

Substituting this into the original equation of motion and diving out the identical exponentials yields

5

Solving for α yields

6

The right hand term is the natural frequency of the system, more generally expressed as a real number:

7

Thus for simplicity the solution can be written as

8

At this point it is not clear which of these two solutions is correct, so let us write the general solution as

9

Because of the complex exponential definition of sines and cosines, we see the beginning of a solution in simply one or the other, but at this point the coefficients are in the way.

These coefficients are determined from the initial conditions. Let us consider these at t=0:

10

Substituting these into our assumed general solution yields

11

The coefficients then solve to

12

It is noteworthy that the two coefficients are complex conjugates of each other.

Since the general solution is written in exponential form, it makes sense that, if the coefficients are to be removed so we can enable a direct solution to a sine or cosine, they too should be in exponential form. Converting the two coefficients to polar form yields

13

Substituting these coefficients into the general solution, we have

14

Factoring out the radical and recognising that the arctangent is an odd function,

15

The quantity in brackets is the complex exponential definition of the cosine, since the two exponents are negatives of each other. The solution can thus be written as

16

If we define

17

the solution is

18

which can be rewritten in a number of ways.

If the dampening is added, the problem can be solved in the same way, but the algebra is a little more complicated, and we will end up additionally with a real exponential (decay) in the final solution.

This derivation demonstrates the power of complex analysis as applied to differential equations even in a simple way.

More examples of this kind of thing are here and here.

Taking the Last Voyage with Newton and Pascal

Don Warrington Fluid Mechanics, Solid Mechanics 1 Comment

He’s not widely known outside of the fields he specialised in, but Adhémar Jean Claude Barré de Saint-Venant (1797-1886, usually known in the Anglophone world as simply Saint-Venant) was one of the premier scientists, engineers and mathematicians of the nineteenth century.  His accomplishments were many and include the following:

  • saint-venant-torsionSuccessful derivation of the Navier-Stokes Equations for a viscous flow before Stokes; these equations are the basis of computational fluid dynamics and the analysis of things that fly.
  • Systematisation and development of methods in the theory of elasticity of solids, including his semi-inverse methods for torsion, important in things such as automobile crank shafts.
  • Methods for the analysis of wave mechanics in bars, which we see in many places, from musical instruments to driven foundation piles.

Saint-Venant was born into a royalist, aristocratic, traditionally Roman Catholic family at a time when it was not safe to be any of these: the French Revolution, at that point stumbling from the Reign of Terror to control of France–and most of Europe–by Napoleon Bonaparte.  It was about the latter where Saint-Venant made a statement about himself that got him into trouble with the “new” Europe.  As described in S. Timoshenko’s History of Strength of Materials:

The political events of 1814 had a great effect on Saint-Venant’s career.  In March of that year, the armies of the allies were approaching Paris and the students of the École Polytechnique were mobilized.  On March 30, 1814, they were moving their guns to the Paris fortification when Saint-Venant, who was the first sargeant of the detachment, stepped out from the ranks with the exclamation: “My conscience forbids me to fight for an usurper…” His schoolmates resented that action very much and Saint-Venant was proclaimed a deserter and never allowed to resume his study at the École Polytechnique.

Saint-Venant’s statement of conscience was at once a political and religious statement, and “progressives” of his day didn’t miss either.  The French, then and now innocent of anti-discrimination legislation or sentiments, made his life miserable. The École Polytechnique was and is France’s premier technical institute of higher learning; getting kicked out of it was the equivalent of, say, being expelled from Princeton or MIT.  He worked in the powder industry for nine years, then was admitted to the École des Ponts et Chausées, where his fellow students shunned him.  He graduated first in his class anyway and began his illustrious career in technical things both theoretical and practical.

In spite of his difficulties within France, his reputation outside of her was another matter.  When François Napoleon Moigno wrote his book on statics, he discovered the following:

He (Moigno) wanted the portion on the statics of elastic bodies to be written by an expert in the theory of elasticity, but every time he asked for the collaboration of an English or a German scientist, he was given the same answer: “You have there, close to you, the authority par excellence, M. de Saint-Venant, consult him, listen to him, follow him.” One of them, M. Ettingshausen, added: “Your Academy of Sciences makes a mistake, a great mistake when it does not open its doors to a mathematician who is so highly placed in the opinion of the most competent judges.” In conclusion Moigno observes: “Fatally belittled in France of which he is the purest mathematical glory, M. de Saint-Venant enjoys a reputation in foreign countries which we dare to call grandiose.”

The French finally broke down and admitted Saint-Venant into the Academy of Sciences in 1868.  He continued his work, much of it from his home, up until the time of his death.  When the President of the Academy announced that passing, he made the following statement:

Old age was kind to our great colleague.  He died, advanced in years, without infirmities, occupied up to the last hour with problems which were dear to him and supported in the great passage by the hopes which had supported Pascal and Newton.

Europeans of the time would not have missed the import of the last statement: Pascal and Newton were Christians, and Saint-Venant was being identified with them as one also.  It was also a statement that Saint-Venant, for all of his achievements and interests which have enriched the world, also had an eternal goal as well.

There’s no evidence that Saint-Venant was ostentatious in his faith walk; descriptions of his life show the contrary.  And–shock to today’s atheist–there’s no evidence that it ever impeded the progress of his research or his thought.  As the statistician and eugenicist Karl Pearson, no friend of Christianity, noted:

The more I studied Saint-Venant’s work, the more new directions it seemed to me to open up for original investigation of the most valuable kind. It suggested innumerable unsolved problems in atomic physics, in impact, in plasticity and in a variety of other branches of elasticity, which do not seem beyond solution, and the solution of which if obtained would be of extreme importance. I felt convinced that a study of Saint-Venant’s researches would be a most valuable directive to the several young scientists, whose recent memoirs shew their interest in elasticity as well as their mathematical capacity. Many of the problems raised by Saint-Venant’s suggestive memoirs were quite beyond my powers of analysis, and I recognised that the most useful task I could undertake, was by a careful account of the memoirs themselves to lead the more competent on to their solution.

saint-venant-wave-equation

The biggest impediment he had to face was the blowback from his stand at the École Polytechnique, and that came from his secularist colleagues.  But, when the end came, all of his colleagues knew where he stood, in this life and the next one.

I spend a lot of time on this site and others talking about sea (and sometimes air) voyages.  And I’ve spent most of my career (and all the academic part of same) in the applied sciences.  But when I take my last voyage into eternity, I want to do it in the same hope of Newton and Pascal–and Saint-Venant and Euler for that matter–namely that which comes from following Jesus Christ out of the grave and into eternal life.

And you should too.

Note: I have a more complete account of Saint Venant’s story with a list of references here.

Determining the Characteristic Polynomial of the Companion Matrix by Use of a Large Matrix

Don Warrington Mathematics Leave a comment

Most proofs of the characteristic polynomial of the companion matrix–an important specific case–proceed by induction, and start with a 2\times2 matrix.  It strikes me that an inductive proof has more force (or at least makes more sense) if a larger matrix is used.  In this case we will use a “large” (numerical analysts will laugh at this characterisation) 10\times10 matrix.

Let us begin by making a notation change. Consider the general polynomial

screenshot_20180414_162144

For this to be monic (one of the requirements for the polynomial in question) we should divide by the last coefficient, thus

screenshot_20180414_162202

Our object is thus to prove that this (or a variation of this, as we will see) is the characteristic polynomial of

screenshot_20180414_162217

The characteristic polynomial of this is the determinant of the following:

screenshot_20180414_162231

(For another application of the characteristic polynomial and the companion matrix, click here.)

To find the determinant, we expand along the first row. But then we discover that only two minors that matter: the one in the upper left corner and the one in the upper right. Breaking this up into minors and cofactors yields the following:

screenshot_20180414_162247

The second matrix, however, is an upper triangular matrix with ones for all of its diagonal entries. Its determinant, therefore, is unity. Also rewriting the coefficient of the second term, we have

screenshot_20180414_162300

or

screenshot_20180414_162311

Repeating this process for the next set of minors and cofactors yields

screenshot_20180414_162322

Note carefully the inclusion of -\lambda in the second term. We can also write this as

screenshot_20180414_162333

Repeating this process until the end, it is easy to see that

screenshot_20180414_162343

or more generally

screenshot_20180414_162357

where n is the degree of the polynomial (and the size of the companion matrix.) If we drop the terms we used to make the polynomial monic, we have at last

screenshot_20180414_162409