Runge–Kutta in Haskell

Last time I implemented a Runge-Kutta method was in December 2001. Back then, I implemented the classical 4th order Runge-Kutta (RK4) method in C. However, in the past decade, whenever I needed a numerical solution of an ordinary differential equation (ODE), I used MATLAB. Now, I would like to use Haskell instead.

Consider the following initial value problem (IVP)

,    

where the vector field and the initial condition are given. Solving an instance of the IVP consists of finding a function that satisfies both the differential equation and the initial condition.

Let be the time-step, and let be the -th time instant, where . Suppose that we solve a given instance of the IVP and obtain a closed-form solution ; we then discretize the closed-form solution to obtain an infinite sequence where . However, not all instances of the IVP have a closed-form solution (or even an analytic solution), and even in the cases where a closed-form solution does exist, it can be quite hard to find such a solution. Therefore, instead of discretizing the closed-form solution, one is tempted to discretize the ODE itself to obtain an approximation of the sequence , which we will denote by . We will call the approximate sequence a numerical solution.

Consider now the following discretized initial value problem (DIVP)

,    

where depends on what particular numerical method one wants to use. We will use the (classical) 4th order Runge-Kutta (RK4) method [1], which is given by

where the variables are given by

In a nutshell, one hopes that by making the time-step “small enough”, the numerical solution  will be “close enough” to , the discretization of the actual solution . This hoping involves a fair amount of faith. Especially so if floating-point arithmetic is used, as underflow, overflow, and round-off error are not to be dismissed.

__________

Haskell implementation of RK4

The following Haskell script implements the RK4 method:

-- define 4th order Runge-Kutta map (RK4) 
rk4 :: Floating a => (a -> a) -> a -> a -> a
rk4 f h x = x + (1/6) * (k1 + 2*k2 + 2*k3 + k4)
            where k1 = h * f (x)
                  k2 = h * f (x + 0.5*k1)
                  k3 = h * f (x + 0.5*k2)
                  k4 = h * f (x + k3)

Let us now test this rather succinct implementation.

__________

Example

Consider the following IVP where , , and

,     .

The closed-form solution of this IVP is . We load the script above into GHCi and compute the numerical solution of the IVP using the higher-order function iterate:

*Main> -- discretize closed-form solution
*Main> let xs = [ exp (-t) | t <- [0,0.1..]]
*Main> -- define function f
*Main> let f x = -x
*Main> -- define time-step h
*Main> let h = 1 / 100000
*Main> -- define initial condition
*Main> let x0 = 1.0
*Main> -- compute numerical solution
*Main> let xs_tilde = iterate (rk4 f h) x0
*Main> -- decimate numerical solution
*Main> let xs_dec = [ xs_tilde !! k | k <- [0,1..], rem k 10000 == 0]

Note that the time-step is . Hence, the numerical solution will have samples in between and , which probably is too fine a discretization. Thus, we decimate the numerical solution by a factor of , so that we are left with only eleven samples in between and , with a uniform spacing of between consecutive samples. We can now compare the discretization of the closed-form solution with the heavily decimated numerical solution, as follows:

*Main> -- print discretized closed-form solution
*Main> take 11 xs
[1.0,0.9048374180359595,0.8187307530779818,
0.7408182206817179,0.6703200460356392,
0.6065306597126333,0.5488116360940264,
0.49658530379140947,0.44932896411722156,
0.4065696597405991,0.36787944117144233]
*Main> -- print decimated numerical solution
*Main> take 11 xs_dec
[1.0,0.9048374180359563,0.8187307530779842,
0.7408182206817214,0.6703200460356442,
0.6065306597126376,0.5488116360940266,
0.49658530379141175,0.44932896411722334,
0.40656965974059905,0.367879441171439]
*Main> -- compute difference of two solutions
*Main> let es = zipWith (-) xs_dec xs
*Main> -- print difference
*Main> take 11 es
[0.0,-3.219646771412954e-15,2.3314683517128287e-15,
3.552713678800501e-15,4.9960036108132044e-15,
4.3298697960381105e-15,2.220446049250313e-16,
2.275957200481571e-15,1.7763568394002505e-15,
-5.551115123125783e-17,-3.3306690738754696e-15]
*Main> -- compute energy of the difference
*Main> sum $ take 11 $ map (^2) es
9.161263559461204e-29

Note how small the values of difference sequence are, of the order of . The implementation of the RK4 method seems to be working. Who needs MATLAB now? ;-)

__________

References

[1] Richard W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 1973.

Life pledged for an ideal

[ Ernst Jünger's famous M1916 Stahlhelm (source) ]

Here is a delirious (and slightly disturbing) passage from the last pages of Basil Creighton’s 1929 translation [1] of Ernst Jünger‘s gut-wrenching masterpiece, In Stahlgewittern [2], perhaps the most beautiful war book ever written:

Now I looked back: four years of development in the midst of a generation predestined to death, spent in caves, smoke-filled trenches, and shell-illumined wastes; years enlivened only by the pleasures of a mercenary, and nights of guard after guard in an endless perspective; in short, a monotonous calendar full of hardships and privation, divided by the red-letter days of battles. And almost without any thought of mine, the idea of the Fatherland had been distilled from all these afflictions in a clearer and brighter essence. That was the final winnings in a game on which so often all had been staked: the nation was no longer for me an empty thought veiled in symbols; and how could it have been otherwise when I had seen so many die for its sake, and been schooled myself to stake my life for its credit every minute, day and night, without a thought? And so, strange as it may sound, I learned from this very four years’ schooling in force and in all the fantastic extravagance of material warfare that life has no depth of meaning except when it is pledged for an ideal, and that there are ideals in comparison with which the life of an individual and even of a people has no weight. And though the aim for which I fought as an individual, as an atom in the whole body of the army, was not to be achieved, though material force cast us, apparently, to the earth, yet we learned once and for all to stand for a cause and if necessary to fall as befitted men.

Hardened as scarcely another generation ever was in fire and flame, we could go into life as though from the anvil; into friendship, love, politics, professions, into all that destiny had in store. It is not every generation that is so favoured.

And if it be objected that we belong to a time of crude force our answer is: We stood with our feet in mud and blood, yet our faces were turned to things of exalted worth. And not one of that countless number who fell in our attacks fell for nothing. Each one fulfilled his own resolve.

__________

Quite a powerful and passionate celebration of the old Preußische Tugenden! I would dare to claim that such extreme self-denial will appear somewhat suicidal and deranged to anyone who is deeply immersed in today’s Zeitgeist.

Interestingly, I can find the text corresponding to the passage above neither in the 1922 German edition [2], nor in Michael Hofmann’s translation [3]. Did Basil Creighton manufacture prose out of thin air? Or is Creighton’s 1929 translation based on a cruder, more brutal, and more nationalistic German edition of In Stahlgewittern? I guess the latter is the case. In any case, if any of you can find the original German text on which Creighton based his 1929 translation, I would be most thankful.

__________

References

[1] Ernst Jünger, Basil Creighton (translator), The storm of steel: from the diary of a German storm-troop officer on the western front, H. Fertig, 1975 (Reprint of the Chatto & Windus, 1929 edition).

[2] Ernst Jünger, In Stahlgewittern, Berlin 1922.

[3] Ernst Jünger, Michael Hofmann (translator), Storm of Steel, Penguin Books, 2004.

Schopenhauer on intellectual vanity

Arthur Schopenhauer on the perils of exhibiting intelligence [1]:

Was für ein Neuling ist doch Der, welcher wähnt, Geist und Verstand zu zeigen wäre ein Mittel, sich in Gesellschaft beliebt zu machen! Vielmehr erregen sie, bei der unberechenbar überwiegenden Mehrzahl, einen Haß und Groll, der um so bitterer ist, als der ihn Fühlende die Ursache desselben anzuklagen nicht berechtigt ist, ja, sie vor sich selbst verhehlet.

You can take a look at the original in [1]. I find the Gothic font rather hard to read and, therefore, I cannot guarantee that I made no mistakes while copying.

__________

Here is T. Bailey Saunders’ translation [2]:

A man must be still a greenhorn in the ways of the world, if he imagines that he can make himself popular in society by exhibiting intelligence and discernment. With the immense majority of people, such qualities excite hatred and resentment, which are rendered all the harder to bear by the fact that people are obliged to suppress–even from themselves–the real reason of their anger.

__________

It would most certainly be instructive to post a longer passage from Saunders’ translation in [3]:

A man must be still a greenhorn in the ways of the world, if he imagines that he can make himself popular in society by exhibiting intelligence and discernment. With the immense majority of people, such qualities excite hatred and resentment, which are rendered all the harder to bear by the fact that people are obliged to suppress—even from themselves—the real reason of their anger.

What actually takes place is this. A man feels and perceives that the person with whom he is conversing is intellectually very much his superior.

He thereupon secretly and half unconsciously concludes that his interlocutor must form a proportionately low and limited estimate of his abilities. That is a method of reasoning—an enthymeme—which rouses the bitterest feelings of sullen and rancorous hatred. And so Gracian is quite right in saying that the only way to win affection from people is to show the most animal-like simplicity of demeanor—para ser bien quisto, el unico medio vestirse la piel del mas simple de los brutos.

To show your intelligence and discernment is only an indirect way of reproaching other people for being dull and incapable. And besides, it is natural for a vulgar man to be violently agitated by the sight of opposition in any form; and in this case envy comes in as the secret cause of his hostility. For it is a matter of daily observation that people take the greatest pleasure in that which satisfies their vanity; and vanity cannot be satisfied without comparison with others. Now, there is nothing of which a man is prouder than of intellectual ability, for it is this that gives him his commanding place in the animal world. It is an exceedingly rash thing to let any one see that you are decidedly superior to him in this respect, and to let other people see it too; because he will then thirst for vengeance, and generally look about for an opportunity of taking it by means of insult, because this is to pass from the sphere of intellect to that of will—and there, all are on an equal footing as regards the feeling of hostility. Hence, while rank and riches may always reckon upon deferential treatment in society, that is something which intellectual ability can never expect; to be ignored is the greatest favor shown to it; and if people notice it at all, it is because they regard it as a piece of impertinence, or else as something to which its possessor has no legitimate right, and upon which he dares to pride himself; and in retaliation and revenge for his conduct, people secretly try and humiliate him in some other way; and if they wait to do this, it is only for a fitting opportunity. A man may be as humble as possible in his demeanor, and yet hardly ever get people to overlook his crime in standing intellectually above them. In the Garden of Roses, Sadi makes the remark:—You should know that foolish people are a hundredfold more averse to meeting the wise than the wise are indisposed for the company of the foolish.

I wonder if “Sadi” refers to the medieval Persian poet Saadi, and if “Garden of Roses” refers to Gulistan.

__________

References

[1] Arthur Schopenhauer, Parerga und Paralipomena: kleine philosophische Schriften, Berlin, 1851.

[2] Arthur Schopenhauer, T. Bailey Saunders (translator), Counsels and Maxims by Arthur Schopenhauer, The Echo Library, 2006.

[3] Arthur Schopenhauer, T. Bailey Saunders (translator), Counsels and Maxims from the Essays of Arthur Schopenhauer, Project Gutenberg, 2004.

Deciding the convexity of semialgebraic sets II

Let denote the set of real symmetric matrices. We introduce a matrix-valued function defined by

The elliptope is defined by

where means that matrix is positive semidefinite. Since the elliptope is defined by the linear matrix inequality (LMI) , we can conclude that  is a spectrahedron and, therefore, we have that is convex.

The elliptope (also known as “inflated tetrahedron”) is the yellow part (surface and interior) of the Cayley’s cubic surface plotted below

[ source ]

__________

is a semialgebraic set

Saying that is positive semidefinite is equivalent to saying that the principal minors of are nonnegative. If , then the three 1st order principal minors (which are the diagonal entries of ) lead to the following three (admittedly uninteresting) inequalities

,     ,    

which we can discard. The three 2nd order principal minors (determinants of submatrices of ) and the 3rd order principal minor (the determinant of ) can be computed using the following Python 2.5 + SymPy script:

from sympy import *

# symbolic variables
x1 = Symbol('x1')
x2 = Symbol('x2')
x3 = Symbol('x3')

# build matrix A(x)
A = Matrix([[1,x1,x2],
            [x1,1,x3],
            [x2,x3,1]])

print "\nMatrix A(x):"
print A

print "\n2nd order principal minors:"
print A.extract([0,1],[0,1]).det()
print A.extract([0,2],[0,2]).det()
print A.extract([1,2],[1,2]).det()

print "\n3rd order principal minor:"
print A.det()

which produces the following output:

Matrix A(x):
[ 1, x1, x2]
[x1,  1, x3]
[x2, x3,  1]

2nd order principal minors:
-x1**2 + 1
-x2**2 + 1
-x3**2 + 1

3rd order principal minor:
-x1**2 + 2*x1*x2*x3 - x2**2 - x3**2 + 1

Hence, the 2nd order principal minors yield the following inequalities

,     ,    

which define the -dimensional cube . The 3rd order principal minor (i.e., the determinant) yields the inequality

.

We can thus conclude that is a basic closed semialgebraic set.

__________

Deciding the convexity of

In my last post, I used quantifier elimination to decide the convexity of two semialgebraic sets in . Let us now decide the convexity of , which is of higher “complexity” than the sets I considered last time. What is the goal? The goal is to find out whether REDLOG can decide the convexity of in less than 10 minutes.

Based on the scripts in my previous post, we write the following REDUCE + REDLOG script to decide the convexity of :

% decide the convexity of the elliptope E3

load_package redlog;
rlset ofsf;

% define predicates in antecedent
p1 := (1-x1^2>=0) and (1-x2^2>=0) and (1-x3^2>=0) and 
      (1-x1^2-x2^2-x3^2+2*x1*x2*x3>=0);
p2 := (1-y1^2>=0) and (1-y2^2>=0) and (1-y3^2>=0) and 
      (1-y1^2-y2^2-y3^2+2*y1*y2*y3>=0);
p3 := (theta>=0) and (theta<=1);

% define predicate in the consequent
z1 := theta*x1+(1-theta)*y1;
z2 := theta*x2+(1-theta)*y2;
q  := (1-z1^2>=0) and (1-z2^2>=0) and (1-z3^2>=0) and 
      (1-z1^2-z2^2-z3^2+2*z1*z2*z3>=0);

% define universal formula
phi := all({x1,x2,x3,y1,y2,y3,theta}, (p1 and p2 and p3) impl q);

% perform quantifier elimination
rlqe phi;

end;

After 30 minutes (!!!) waiting for the decision procedure to terminate, I killed the reduce.exe process (which was eating all my machine’s CPU cycles and RAM memory). This is disappointing…

Deciding the convexity of semialgebraic sets

Given a basic closed semialgebraic set [1], , defined by

where and , how can one decide if set is convex? Can one use quantifier elimination?

__________

Convexity

Let us recall the definition of convexity [2]:

Definition: A set is convex if the line segment between any two points in lies in , i.e., if for any and any with , we have .

Let us introduce a predicate , defined by

so that we can write in the more parsimonious form

.

Hence, writing is equivalent to asserting that . From the definition above, it follows that set is convex if and only if the following universally quantified formula

where range over and ranges over , evaluates to . Since are polynomials, the formula above can be decided using quantifier elimination software like QEPCAD or REDLOG. Unfortunately, decidability does not imply that real quantifier elimination will decide the formula in an expedite manner. Let us now consider two instances of the decision problem under study.

__________

Example #1

Consider the following semialgebraic set

which we depict below

Visual inspection of the plot above allows us to conclude that set is convex. However, relying on the human visual system is extremely limiting. Therefore, we would like to automate the decision. We decide the convexity of via quantifier elimination using the following REDUCE + REDLOG script:

% decide the convexity of a semialgebraic set

load_package redlog;
rlset ofsf;

% define predicates in antecedent
p1 := (x1-x2^2+3>=0) and (-x1-x2^2+3>=0);
p2 := (y1-y2^2+3>=0) and (-y1-y2^2+3>=0);
p3 := (theta>=0) and (theta<=1);

% define predicate in the consequent
z1 := theta*x1+(1-theta)*y1;
z2 := theta*x2+(1-theta)*y2;
q  := (z1-z2^2+3>=0) and (-z1-z2^2+3>=0); 

% define universal formula
phi := all({x1,x2,y1,y2,theta}, (p1 and p2 and p3) impl q);

% perform quantifier elimination
rlqe phi;

end;

This script returns the truth value and, hence, is, indeed, convex. However, it took REDLOG approximately 60 seconds to decide the convexity of , which is considerably longer than any of my previous experiments with REDLOG. Since we are performing quantifier elimination on a formula with five universal quantifiers, I am not incredibly surprised that it takes a while to obtain an answer.

__________

Example #2

Consider now the following semialgebraic set

part of which we depict below

Visual inspection of the plot above allows us to conclude that set is not convex. We decide the convexity of using the following REDUCE + REDLOG script:

% decide the convexity of a semialgebraic set

load_package redlog;
rlset ofsf;

% define predicates in antecedent
p1 := (x1-x2^2+3>=0) and (x1-x2^2+1<=0);
p2 := (y1-y2^2+3>=0) and (y1-y2^2+1<=0);
p3 := (theta>=0) and (theta<=1);

% define predicate in the consequent
z1 := theta*x1+(1-theta)*y1;
z2 := theta*x2+(1-theta)*y2;
q  := (z1-z2^2+3>=0) and (z1-z2^2+1<=0); 

% define universal formula
phi := all({x1,x2,y1,y2,theta}, (p1 and p2 and p3) impl q);

% perform quantifier elimination
rlqe phi;

end;

This script returns the truth value and, hence, is, indeed, not convex. Interestingly, it took REDLOG less than one second to decide the formula. That is over two orders of magnitude faster than in example #1.

__________

References

[1] Markus Schweighofer, Describing convex semialgebraic sets by linear matrix inequalities, International Symposium on Symbolic and Algebraic Computation, Korea Institute for Advanced Study, Seoul, July 2009.

[2] Stephen Boyd, Lieven Vandenberghe, Convex Optimization.