Earth Cubed

Distributed Climate Science and Computing

Laplace Transform Via Limits

Quite a while back I had a brief interest in wheather the Laplace transform made numeric sense. I was interested in this for its own sake:

and because if we could numericaly understand the Laplace transform it might tell us more about a dynamic system then a fourier transform.

Memory of my failures in obtaining the desired results led me to believe that the Laplace transform was not Riemann integratable. However when chalanged on this in another forum:

I rethought a suggested solution and decided that I was perhaps partly incorrect. We can derive the integral from sums & limits bit it converges in a very strange way. The sum converges exactly where the geometric series diverges. All that is left prior to cancellation is a ratio of infinitesimally small quantities.

The summation suggested by andrewk for the laplace transform was:

(1) – F(s)= \int_0^{oo} f(t) dt = lim_{m->oo} lim_{n->oo} \sum_{k=1}^{n} f(km/n) \frac{m}{n}

It looks simple enough that I may have or should have considered it but it is not entirely obvious for how the double sum should work. Moreover, wikipedia discusses some problems with trying to use Riemann integration with improper integrals.

Needless to I am able to use a Riemann like result and get the right answer. Whether the approach is mathematically sound is an entirely different question. In the above equation if we are trying to integrate an exponentail function of the form

(2) – f(x)=exp(- \alpha t)

(3) – L( f(x)) = lim_{m->oo} lim_{n->oo} \sum_{k=1}^{n} e^{-skm/n} exp(- \alpha km/n) (m/n) 

combining terms

(4) – L[f(x)] =lim_{m->oo}lim_{n->oo}\sum_{k-1}^n exp(-(s+alpha)km/n) (m/n) 

Expressing as a geometric series:

(5) – L[f(x)]=lim_{m->oo}lim_{n->oo}\sum_{k-1}^n \left( e^{-(s+alpha)m/n} \right)^k (m/n) 

Using the geometric series:

(6) – L[f(x)]=lim_{m->oo}lim_{n->oo} \frac{1- \left( e^{-(s+alpha)m/n \right)^n}}{1-e^{-(s+alpha)m/n}} \frac{m}{n} 

Now here is where the Voodoo comes in we need to let m/n approach zero and m approach infinity. Rather them attempting to do this formally lets just cut to the chase.

when m/n approaches zero we can we can approximate:$latex e^{-(s+alpha)m/n} $ as 1 - (s+alpha) (m/n) 

and as m approach infinity \left( e^{-(s+alpha)m/n \right}^n   approaches zero.

Making these substitution into equation (6) we get:

(7) – L[f(x)]=lim_{m->oo}lim_{n->oo} \frac{1}{1-(1+(s+alpha)m/n)} \frac{m}{n}

which further simplifies to

(8) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha)m/n)} x m/n

cancelling m/n from the numerator and denominator (

(9) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha))}

and now the limits are irrelivant:

(10) – L{f(x)}=lim_{m->oo}lim_{n->oo} \frac{1}{(s+alpha))}


July 14, 2012 Posted by | Uncategorized | Leave a comment

log(CO2) and Scary Graphs

After reading Aurther’s blog post “New Congressional Budget Office Report on Climate Change” I got curious as to what the temperature response to CO2 would look like if it is actually logarithmic.

I descovered (what I should have realized from the start) that the curve is concave up but will converge to a linear curve for large T

It is commonly believed that the response of the earth to greenhouse gasses is logarithmic. I heard people suggest on other forums that this view was obtained empirically though climate models and perhaps more specifically radiative transfer models.

For analytic justification one would derive an expression for how the spectral width of an absorption peak grows with CO2 concentration (as I’ve done here). However, I am not convinced that this is sufficient justification as there will always be some radiative transfer at a given wave length even if the majority of it is absorbed over a very short distance. This is because the temperature gradient produced by the lapse rate produces an outward radiative flux that would exceed radiative feedback. I discuss this in more detail in my post:
Tropospheric Feedback

The logarithmic response is important because it is a type of saturation, in other words the more CO2 that is added to the atmosphere the less effective the next unit of CO2 will be in contributing to the warming. What I learned from reading Arther’s blog is that the current CO2 levels have not yet overwhelmed the natural levels of CO2.  This can be seen in the following graph:

More specificity from about 1000-1800 the CO2 concentration in the atmosphere stayed around 280 ppm. The following graph is more useful for measuring the current growth in CO2 concentration:

This graph is surprisingly very linear. If the growth in CO2 is truely exponential then it should be possible to estimate in from the slope on this graph which is given as 1.4203 PPM per year. For an exponential function:

y=y_o+Aexp(\lambda (t-t_o))
The derivative is:
y'=A\lambda exp( \lambda (t-t_o))
And the second derivative is:
y''=A\lambda^2 exp( \lambda (t-t_o))
The second derivative was taken because two equations are needed to find both A \lambda can be found.

The site also where I obtained the above figures gives a quadratic fit which can be used to estimate the first and second derivatives:


Therefore at year 2007 the first derivative is given by:


and the second derivative is 0.0119942


A exp( \lambda (t-t_o)) \lambda =1.9961
A exp( \lambda (t-t_o)) \lambda^2 =0.0119942

Dividing the second equation by the first:


From this the doubling time can be obtained as follows:

ln(2)/\lambda=ln(2)/0.0060=115.5245 years.

y_o is taken to be the base level of CO2. That is:


A third equation can now be found as:




If t_o is taken to be 1800 this gives:


which suggests the CO2 growth rate has decreased over the last 200 years.

The CO2 is estimated to follow this function:

CO_2(t)=280+31.8931 exp(0.006 (t-1800))

The question now is how does this growth rate in CO2 effect the temperature. There are several estimates for the sensitivity of the climate to changes in CO2. Lucia’s one box model “lumpy” suggest a sensitivity of:

1.7 degrees Celsius per CO2 doubling.
The IPCC estimates the lower bound for sesitivity to be:
1.5 degrees Celsius per CO2 doubling (see CO2 Climate Sensitivity)
Isaac M. Held suggests a climate sensitivity of about 2.8C/CO2 doubling. See:

Selected recent papers on climate sensitivity:

Here is what wikipedia has to say:

In Intergovernmental Panel on Climate Change (IPCC) reports, equilibrium climate sensitivity refers to the equilibrium change in global mean near-surface air temperature that would result from a sustained doubling of the atmospheric (equivalent) CO2 concentration. This value is estimated, by the IPCC Fourth Assessment Report (AR4) as likely to be in the range 2 to 4.5°C with a best estimate of about 3°C, and is very unlikely to be less than 1.5°C. Values substantially higher than 4.5°C cannot be excluded, but agreement of models with observations is not as good for those values. This is a slight change from the IPCC Third Assessment Report (TAR), which said it was “likely to be in the range of 1.5 to 4.5°C” [1]. AR3 defined climate sensitivity alternatively in systematic units, equilibrium climate sensitivity refers to the equilibrium change in surface air temperature following a unit change in radiative forcing and is expressed in units of °C/(W/m2) or equivalently K/(W/m2). In practice, the evaluation of the equilibrium climate sensitivity from models requires very long simulations with coupled global climate models, or it may be deduced from observations. Therefore the 2007 AR4 renamed the alternative climate sensitivity to climate sensitivity parameter adding a new definition of effective climate sensitivity which is “a measure of the strengths of the climate feedbacks at a particular time and may vary with forcing history and climate state”.

The logarithmic law of CO2 forcing is given as:

\Delta T=k*log_2(CO_2(t_2)/CO_2(t_1))=\frac{k}{log_n(2)}log_n((CO_2(t_2)/CO_2(t_1))

where k is the CO2 sensitivity for doubling CO2

I plotted this function for several values of the doubling sensitivity k

The labels on the right hand side of the plot are the climate sensitivities for each curve. This is actually a considerably smaller response then one would expect given the doubling time is around 100 years. However, while this is a nearly sufficient time for the exponential part of the curve to double,  the CO2 only increases by a factor of 1.2 since at 1997 the CO2 concentration is 384 ppm and in 2100 the CO2 concentration is projected by this fit to be  473ppm. This reduces the expected response by a factor of:


Notice that if the sensitivity ranges given by the IPCC are used then with this fit to CO2 emission growth is around 0.5-1.2 degrees which is hardly the doomsday scenario shown in the following graph which was posted on Aurther’s blog.

As a final note, the MATLAB code I used to produce the above graph is:

clear all
CO2= @(t)280+31.8931*exp(0.006*(t-1800))
K=[1.5 4 7 11 15 21 25 29 33]
for i=1:length(K)
  AXIS([1997 2100 0 10])
  hold on;
ylabel('Temperature Change in Degrees Celcius')
hold off,

September 19, 2009 Posted by | Uncategorized | 1 Comment

Numeric Solutions to The Heat Equation

I have been reading a lot on Lucia’s blog about two box models which are essentially an approximation of the heat equation with basis functions which are constant over a box.

The heat equation is given by:

\frac{\partial u}{\partial t} -k\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\right)=0

or equivalently:

\frac{\partial u}{\partial t} = k \nabla^2 u

The fundamental solutions or Greens functions (also see main body theory) are of the form:

\Phi(\mathbf{x},t) = \Phi(x_1,t)\Phi(x_2,t)\dots\Phi(x_n,t)=\frac{1}{(4\pi k t)^{n/2}}e^{-\mathbf{x}\cdot\mathbf{x}/4kt}.

This suggests my choice of a negative exponential basis in my post (Lagrangian Mechanics and the Heat Equation) was not two bad a choice although, Guasian functions will decay slightly faster then negative exponentials.

Not all solutions are based on on fundamental solutions for instance in the post (Approximations used in Crank-Nicolson method for solving PDEs numerically) I read that the Crank-Nicolson method was the standard method of soliving the Heat equation numericaly.

For instance in 1-D

\frac{\partial u}{\partial t}=k \frac{\partial^2u}{\partial x^2}

the Cranck Nicholson Method is given by:

\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = \frac{k}{2 (\Delta x)^2}\left((u_{i + 1}^{n + 1} - 2 u_{i}^{n + 1} + u_{i - 1}^{n + 1}) + (u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n})\right)

It should be noted that this method produces a difference equation. The values at the next time step can be solved for analytically, using cramer’s rule (see also Invertible matrix, analytic solutions). The frequency domain characteristics can be explored using the z transform. Where the frequency response is given by evaluating the z transform along the unit circle.

Also note that numeric error can be reduced when computing future time steps by either recursive squaring:
or by using matrix decomposition.

For other numeric method of solve partial differential equations see (Numerical partial differential equations), which I posted in the thread (Preparation for PDEs).

Further with regards to crank Nicholson, there is no time dependency on the right hand side of the equations so other methods, can be used to descretize the heat equations such as using Laplace transforms or the matrix exponential.

With regards to Lucia’s blog

My understanding as posted in (Arthur’s Case 2 (I think)) that the main focus of Lucia’s blog posts is to test the model chosen by Tamino:

lucia (Comment#19822) September 12th, 2009 at 9:23 pm

What I mean is– when testing the two box model, you don’t switch to the diffusive model even if it’s more inherently sensible and intelligent. That’s because to test “X” you must test “X”. You can’t test “Y” even if “Y” seems more likely to pass the test.

This is fine but I think that a wider discussion is warranted about how this model is just a simplified version of the heat equations and what principles of modeling and differential equations can be useful to obtain better solutions.

September 15, 2009 Posted by | Math | 3 Comments

Coriolis Forces in Hopkins and Simmons Vorticity Equation

In the thread vector operations in Hopkins and Simmons, I compute the components of the curl as:

\left[ \nabla \times \vec A \right]_{\sigma}= \frac{1}{r\mu}\frac{\partial A_{\lambda}}{\partial \mu}+ \frac{RT}{\sigma m g}\frac{\partial A_{\mu}}{\partial \lambda}=\left[ \nabla \times \vec A \right]_{\sigma,\mu}+\left[ \nabla \times \vec A \right]_{\sigma,\lambda}

\left[ \nabla \times \vec A \right]_{\mu}=  -\frac{RT}{\sigma m g}\frac{\partial A_{\sigma}}{\partial \lambda}- \frac{\sqrt{1-\mu^2}}{r} \frac{\partial A_{\lambda}}{\partial \sigma}=\left[ \nabla \times \vec A \right]_{\mu,\lambda}+\left[ \nabla \times \vec A \right]_{\mu,\sigma}

\left[ \nabla \times \vec A \right]_{\lambda}=  \frac{\sqrt{1-\mu^2}RT}{r}\frac{\partial A_{\mu}}{\partial \sigma}-\frac{RT}{r\mu} \frac{\partial A_{\sigma}}{\partial \mu}=\left[ \nabla \times \vec A \right]_{\lambda,\sigma}+\left[ \nabla \times \vec A \right]_{\lambda,\mu}

In my post Coriolis forces in Hopkins and Simmons I compute the coriolis force as:

\boldsymbol{\Omega \times v} =  \begin{pmatrix} \mu U \\  \pm \sqrt{1-\mu^2} U \\ \mp \sqrt{1-\mu^2} V - \mu  W \end{pmatrix} \  = \begin{pmatrix} A_{\sigma} \\ A_{\mu} \\ A_{\lambda} \end{pmatrix} \

And the partial derivatives are given by:

\frac{\partial A_{\lambda}}{\partial \mu}=\frac{\partial}{\partial \mu}\left(\mp \sqrt{1-\mu^2} V - \mu  W\right)=\pm \frac{2\mu}{\sqrt{1-\mu^2}}V \mp \sqrt{1-\mu^2} \frac{\partial V}{\partial \mu} - W-\frac{\partial W}{\partial \mu}

\frac{\partial A_{\mu}}{\partial \lambda}=\frac{\partial }{\partial \lambda}\left(\pm \sqrt{1-\mu^2} U\right)=\pm \sqrt{1-\mu^2} \frac{\partial U}{\partial \lambda}

\frac{\partial A_{\sigma}}{\partial \lambda}=\frac{\partial }{\partial \lambda}\left( \mu U \right)=\mu \frac{\partial U}{\partial \lambda}

\frac{\partial A_{\lambda}}{\partial \sigma}=\frac{\partial}{\partial \sigma}\left(\mp \sqrt{1-\mu^2} V - \mu  W\right)=\mp \sqrt{1-\mu^2} \frac{\partial V}{\partial \sigma} - \mu  \frac{\partial W}{\partial \sigma}

\frac{\partial A_{\mu}}{\partial \sigma}=\frac{\partial }{\partial \sigma}\left(\pm \sqrt{1-\mu^2} U\right)=\pm \sqrt{1-\mu^2} \frac{\partial U }{\partial \sigma}

\frac{\partial A_{\sigma}}{\partial \mu}=\frac{\partial }{\partial \mu}\left( \mu U \right)=\mu \frac{\partial U}{\partial \mu}

I’ll derive the rest of this later but this doesn’t seem to be the form of the prognostic equation used by Hopkins and Simmons.

Divergence Free Flow

In the post Vector Operations in Hopkins and Simmons I derived the divergence operator as follows:

\nabla \cdot  = \frac{1}{h_{\mu}} \frac{\partial }{\partial \mu} + \frac{1}{h_{\lambda}}\frac{\partial}{\partial \lambda}+\frac{1}{h_{\sigma}}\frac{\partial}{\partial \sigma}=\frac{1-\mu^2}{r} \frac{\partial}{\partial \mu} + \frac{1}{r \mu}\frac{\partial}{\partial \lambda}-\frac{RT}{\sigma mg}\frac{\partial}{\partial \sigma}

If the divergence of the velocity equals zero then:

\nabla \cdot \begin{pmatrix} W \\ V \\ U \end{pmatrix}\ =\frac{1-\mu^2}{r} \frac{\partial V}{\partial \mu} + \frac{1}{r \mu}\frac{\partial U}{\partial \lambda}-\frac{RT}{\sigma mg}\frac{\partial W}{\partial \sigma}=0

Which implies:

\frac{\partial W}{\partial \sigma}=\frac{(1-\mu^2)\sigma mg}{r RT} \frac{\partial V}{\partial \mu} + \frac{\sigma mg}{r \mu \sigma mg}\frac{\partial U}{\partial \lambda}

September 12, 2009 Posted by | GCM (General Circulation Models) | 1 Comment

The Cross Product in Non Orthogonal Coordinate Systems

The form of the cross product I’ve shown in my post Coriolis Forces is:

\boldsymbol{\Omega \times v} =  \begin{pmatrix} \Omega_y v_z - \Omega_z v_y \\ \Omega_z v_x - \Omega_x v_z  \\ \Omega_x v_y - \Omega_y v_x  \end{pmatrix}

The components of this cross product can be written as follows:

\left[ \boldsymbol{\Omega \times v} \right]_x=  \left[ \begin{matrix}\Omega_x&\Omega_y&\Omega_z \end{matrix} \right]\left[ \begin{matrix}0&0&0 \\ 0&0&1 \\ 0&-1&0\end{matrix} \right]\left[ \begin{matrix}v_x\\v_y\\v_z \end{matrix} \right]

\left[ \boldsymbol{\Omega \times v} \right]_y=  \left[ \begin{matrix}\Omega_x&\Omega_y&\Omega_z \end{matrix} \right]\left[ \begin{matrix}0&0&-1 \\ 0&0&0 \\ 1&0&0\end{matrix} \right]\left[ \begin{matrix}v_x\\v_y\\v_z \end{matrix} \right]

\left[ \boldsymbol{\Omega \times v} \right]_z=  \left[ \begin{matrix}\Omega_x&\Omega_y&\Omega_z \end{matrix} \right]\left[ \begin{matrix}0&1&0 \\ -1&0&0 \\ 0&0&0\end{matrix} \right]\left[ \begin{matrix}v_x\\v_y\\v_z \end{matrix} \right]

We will abbreviate these relationships as follows:

\left[ \boldsymbol{\Omega \times v} \right]_x=\Omega^T[\times]_xv

\left[ \boldsymbol{\Omega \times v} \right]_y=\Omega^T[\times]_yv

\left[ \boldsymbol{\Omega \times v} \right]_z=\Omega^T[\times]_zv

Now define the coordinate transform:

\boldsymbol{r'}=T \boldsymbol{r}
\boldsymbol{r}=\left[\begin{matrix}  x \\ y \\ z \end{matrix} \right]

Then the cross product components can be written as follows:

\left[ \boldsymbol{\Omega \times v} \right]_x=(\Omega^TT^T)((T^{-1})^T[\times]_xT^{-1})(Tv)

\left[ \boldsymbol{\Omega \times v} \right]_y=(\Omega^TT^T)((T^{-1})^T[\times]_yT^{-1})(Tv)

\left[ \boldsymbol{\Omega \times v} \right]_z=(\Omega^TT^T)((T^{-1})^T[\times]_zT^{-1})(Tv)

Now Right multiplying the matrix \left[ \boldsymbol{\Omega \times v} \right] by the transform T gives:

\left[ \boldsymbol{\Omega \times v} \right]_{x'}=\sum_{k \in \{x,y,z\}}T_{1,k}(\Omega^TT^T)((T^{-1})^T[\times]_kT^{-1})(Tv)

\left[ \boldsymbol{\Omega \times v} \right]_{y'}=\sum_{k \in \{x,y,z\}}T_{2,k}(\Omega^TT^T)((T^{-1})^T[\times]_kT^{-1})(Tv)

\left[ \boldsymbol{\Omega \times v} \right]_{z'}=\sum_{k \in \{x,y,z\}}T_{3,k}(\Omega^TT^T)((T^{-1})^T[\times]_kT^{-1})(Tv)

Which can be written in this form:
\left[ \boldsymbol{\Omega \times v} \right]_{x'}=\Omega'^T[\times]_{x'}v'

\left[ \boldsymbol{\Omega \times v} \right]_{y'}=\Omega'^T[\times]_{y'}{v'}

\left[ \boldsymbol{\Omega \times v} \right]_{z'}=\Omega^T[\times]_{z'}{v'}


\left[[\times]_{x'}\right]_{i,j}=\sum_{k \in \{x,y,z\}}T_{1,k}[((T^{-1})^T[\times]_kT^{-1})]_{i,j}

\left[[\times]_{y'}\right]_{i,j}=\sum_{k \in \{x,y,z\}}T_{2,k}[((T^{-1})^T[\times]_kT^{-1})]_{i,j}

\left[[\times]_{z'}\right]_{i,j}=\sum_{k \in \{x,y,z\}}T_{3,k}[((T^{-1})^T[\times]_kT^{-1})]_{i,j}


September 7, 2009 Posted by | Math | 2 Comments

Lagrangian Mechanics and The Heat Equation

I have been noticing a lot of discussion on Lucia’s blog about an energy Balance Climate Model.

Two Box Model: Algebra for ver. 1 test.
Two Box Model: Now assuming surface temperature are ‘mixed’ values.
Two Box Model: Rough idea how to obtain parameters.
Two Box Models & The 2nd Law of Thermodynamics

Which started based on a post by Tamino, called “Two Box”

The biggest criticism I’ve scene about these models is that there is very little physics which is used to form these models and consequently they won’t capture much of the nonlinear dynamics.

The basic premises of the model is that the time lags are associated with the heat capacity of the ocean. Consequently, the ocean and some of the atmosphere is partitioned into boxes and the model tries to determine how these boxes are coupled.

The end result though is the time legs end up being associated with eignvectors which, project onto both boxes.  Consequently the modes of the system do not coincide with the boxes.  Moreover, it is desirable to minimize the number of model parameters which need to be fit and therefore principles of physics, should be used to derive a more realistic “Two box model” or should I say more generally “Two mode model” as the modes do not coincide with the boxes.

The most simple equation with regards to the transfer of energy is the heat equation:

1)               \frac{\partial u}{\partial t} -k\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\right)=0

The heat equation is essentially a diffusion equation based on Brownian Motion. But we can adjust the constant K to try and account for other heat exchange processes.

In one  dimension the heat equation can be written as:

2)                 \frac{\partial u}{\partial t} -k\frac{\partial^2u}{\partial y^2}=0

In this equation the time derivative only depends on the spatial derivatives.  Keeping a mind on the dynamics, the higher the frequency that the climate forcing is, the more shallow it will penetrate into the ocean because the more ocean you include the greater the thermal inertia (heat capacity) and consequently the greater the damping.  Thus a good fit for the temperature response at a given frequency as a function of depth, may be an exponential function of the depth.

Thus, from this observation and because the mathematics are fairly simple consider a set of basis functions which the special component is exponentially decreasing from the surface.  That is let:

3)                u=\sum a_i exp(-\lambda_i y)

Plugging this into the heat equation one gets:

4)                \frac{du}{dt}=\sum \frac{d ( \ a_i)}{dt} exp(-\lambda_i )=k \lambda_i^2\sum a_i exp(-\lambda_i y)


The constant k(y) represents how quickly heat diffuses in the ocean at a given depth. Equating terms on the heat equation gives:

5)                k\lambda_i^2a_i-\frac{d (a_i)}{dt}=0

Now let there be a function c(y) , and let the quanity being diffused in the heat equation be the temperature. Then the energy constraint is given by:

6)               E=\int_0^{d} c(y)\sum a_i exp(-\lambda_i y)dy

This constraint makes one of the generalized coordinates a_i redundant.

Differentiating (5) with respect to time, gives:

7)             \frac{dE}{dt}=\int_0^{d} \sum c(y)\frac{d(a_i)}{dt} exp(-\lambda_i y)dy

\frac{dE}{dt}   is determined by the forcing. Now some notes, with regards to the coordinate system:

1)I suggest taking the as the redundant basis function, the one that has the smallest value \lambda .
2)Since this model is a linearization (or approximation), I suggest using as the temperature variable at a given depth, the amount by which it exceeds the average temperature at that depth.  That is model the temperature anomaly instead of the actual temperature.
3)Given two then we are looking at the change in the heat transfer from the mean instead of the actual heat transfer.
4)These are only suggestions.

Now the Lagrangian, is essentially a cost function.  In our cost function, we want to minimize the error in the derivatives (the square of equation (5) integrated over y) and as an additional constraint, I want to minimize the change in entropy with respect to y.

When a parcel of water rises because it is hotter then the surrounding and therefore less dense it will expand with constant entropy (adiabatic expansion) if it exchanges no heat with the surroundings. Therefore, the heat induced ocean currents will act in order to try and minimize the entropy change between adjacent regions. (Note the actual effect of this process is to increase entropy).

From thermodynamics:

8 ) dU=T\,dS-p\,dV


9) dS=(dU-p\,dV)/T=\frac{1}{T}(\frac{1}{C(y)}dT-p(y) \frac{dV}{dP}dP)

from which we can get:

10) \frac{dS}{dy}=\frac{1}{T}(\frac{1}{C(y)}\frac{dT}{dy}-p(y) \frac{dV}{dP}\frac{dP}{dy})

Substituting into (10) the expression for temperature (equation (3)) equation (10) gives:

11) \frac{dS}{dy}=-\frac{1}{T}(\frac{1}{C(y)}\sum a_i\frac{-1}{\lambda_i} exp(\lambda_i y)+p(y) \frac{dV}{dP}\frac{dP}{dy})

The entropy cost function is:

12) \int_o^d (\frac{dS}{dy})^2dy=\int_0^d\frac{1}{T^2}(\frac{1}{C(y)}\sum a_i\frac{1}{\lambda_i} exp(\lambda_i y)+p(y) \frac{dV}{dP}\frac{dP}{dy})^2dy

The Lagrangian is given by:

13) \mathcal L(a_1,...a_n, ...;\dot a_1,..,\dot a_n, ...;t)\,.= \alpha \int_o^d (\frac{dS}{dy})^2dy+\sum(k\lambda_i^2a_i-\frac{d (a_i)}{dt})^2

Where \alpha is a stiffness parameter for the entropy.

Equation 13) is subject to the constraints for some a_j  (These are the energy constraints mentioned above.

14) \frac{d(a_j)}{dt}=\frac{\frac{dE}{dt}-\int_0^{d} \sum_{i \neq j} c(y)\frac{d(a_i)}{dt} exp(-\lambda_i y)dy}{\int_0^{d}c(y)exp(-\lambda_j y)dy}

15) a_j=\frac{E-\int_0^{d} c(y)\sum_{i \neq j} a_i exp(-\lambda_i y)dy}{\int_0^{d} c(y) exp(-\lambda_i y)dy}

Now the Euler Lagrange Equation is used to obtain the dynamics (differential equations) from the Lagrangian

16) \left(\frac{\mathcal{L}(a_1,\ldots,a_n,\ldots,\dot a_1,\ldots,\dot a_n,\ldots;t)}{\partial a_i}\right) - \left(\frac{d}{dy} \frac{\mathcal{L}(a_1,\ldots,a_n, \ldots,;\dot a_1,\ldots,\dot a_n, \ldots;t)}{\partial \dot a_i}\right) = 0 \quad \textrm{for } i = 1, \dots, n, i \neq j

This will give a second order differential equation and it might be necessary to use some linear algebra to rearrange the equation. You can convert this into a first order differential equation by using the Hamiltonian form.  If only three basis  functions are used then one gets “two mode model”. There are no restrictions on the number of bais functions used.

September 1, 2009 Posted by | GCM (General Circulation Models) | 7 Comments

Laplace Transform of f(t) Related to smoothed f(t)?

When reading (Comment#18839) I started to wonder if there was a relationship between the Fourier Transform of a smoothed signal and the Laplace transform.  I assumed there was a relationship (Comment#18854). After further derivation, I recommenced that if the goal is to derive the Laplace transform from the Fourier transom of the filtered signal:

1) The signal be properly windowed.

2) The FFT of the windowed Fourier Transform, needs to be compensated for the frequency effects that resulted from the low pass filter.

Weather it is a good idea to compute the Laplace transform from a windowed FFT of a filtered signal is outside of the scope of this thread (but feel free to comment bellow) .

The Laplace transform is given by:

1)           (\mathcal{L}f)(s) = \int_{0^-}^\infty e^{-st}f(t)\,dt

The Fourier transform is given by:

2)           \hat{f}(\omega) := \int_{-\infty}^{\infty} f(x)\ e^{- \omega i t }\,dt,

The Two Sided Laplace Transform is given by:

3)           \mathcal{B} \left\{f(t)\right\} = F(s) =  \int_{-\infty}^{\infty}  e^{-st} f(t) dt.

Therefor the Fourier transform is the two sided Laplace transform evaluated at s=i\omega

Returning to the one sided Laplace transform:

4)            (\mathcal{L}f)(s) = \int_{a}^b e^{-st}f(t)\,dt= \int_{0^-}^\infty e^{-(\alpha+iw)t}f(t)\,dt

5)            (\mathcal{L}f)(s) = \int_{a}^b e^{-st}f(t)\,dt= \int_{0^-}^\infty e^{-iw t}f(t)e^{- \alpha t}\,dt

Integrating by Parts:


6)            u=e^{- iw t}
7)           dv=f(t)e^{- \alpha t}\,dt

8)            V=\int_{a}^tf(t')e^{-\alpha t'}\,dt'=e^{-\alpha t}\int_{t'=a}^tf(t')e^{\alpha (t-t')}\,dt'

where the low pass filtered version of f(t):

9)            g(t)=f(t)*e^{\alpha t} =e^{-\alpha t} \int_{a}^tf(t')e^{\alpha (t-t')}\,dt'

and is the convolution of f(t) and the impulse response of a filter (or atleast aproximatly so) with bandwidth \alpha .

Plugging this result into integration by parts gives:

10)         (\mathcal{L}f)(s)=[e^{i \omega t}g(t)]_a^b-\int_a^b-i \omega e^{-i \omega t}e^{-\alpha t}g(t)dt

or equivalently:

11)          (\mathcal{L}f)(s)=g(b)e^{-i\omega b}-g(a)e^{-i \omega a}+i\int_a^b \omega e^{-i \omega t}e^{-\alpha t}g(t)dt

The first two terms show how the endpoints chosen effect the transfrom.  These two terms will cancel for a given frequency if the distance between the endpoints is some multiple of the period. The last term is the Fouier transform of the smoothed function with the frequencies weighted by \omega and using a windowing function e^{-\alpha t}
(note the multiple i is there because the Fourier transform variable \omega is the Laplace transform variable but rotated by 90 degrees.)

The effect of the windowing function is to smooth the frequency response. This is because multiplication in the time domain is equivalent to convolution in the frequency domain. The following Fourier transform relationship is useful (relationship 205):

12)           \mathcal{F}(e^{- \alpha t} u(t)) \,=\frac{1}{\sqrt{2 \pi} (\alpha + i \omega)}=\frac{1}{\sqrt{2 \pi (\alpha^2 + \omega^2)}}e^{(-i \ tan^{-1}(\alpha,\omega))}

Note, that if a non causal filter was used for the smoothing the relationship is much simpler.

13)           \mathcal{F}({e}^{-a|t|}) \,=\sqrt{\frac{2}{\pi}} \cdot \frac{a}{a^2 + \omega^2}

In both cases to properly deal with the end points the time shifting property of the Fourier transform is needed (relationship 102):

14)           f(x - t)\, =e^{- i a \omega} \hat{f}(\omega)\,

Applying this property to the last two relationships gives:

15)           \mathcal{F}(e^{- \alpha (t-a)} (u(t-a)-u(t-b)) \,=\frac{e^{-ia \omega}-e^{\alpha a-ib\omega}}{\sqrt{2 \pi} (\alpha + i \omega)}

16)           \mathcal{F}({e}^{-a|(t-b)|}) \,=\sqrt{\frac{2}{\pi}} \cdot \frac{a}{a^2 + \omega^2}e^{-i\omega\frac{1}{2}(a+b)}

Strictly dealing with the case where a causal filter is used and applying the rule for the Fourier transform of a convolution (Rule 109) we obtain:

17)           (\mathcal{L}f)(s)=g(b)e^{-i\omega b}-g(a)e^{-i \omega a}+i\omega\frac{1}{\sqrt{2\pi}}\mathcal{F}(g(t))*\left( \frac{e^{-ia \omega}-e^{\alpha a-ib\omega}}{\sqrt{2 \pi} (\alpha + i \omega)} \right)

of equivalently:

18)           \mathcal{L}(f(y))=g(b)e^{-i\omega b}-g(a)e^{-i \omega a}+i\omega\frac{1}{\sqrt{2\pi}}\mathcal{F}(g(t))*\left( \frac{e^{-ia \omega}-e^{\alpha a-ib\omega}}{\sqrt{2 \pi(\alpha^2 +\omega^2)}}e^{ -i tan^{-1}(\alpha,\omega) } \right)

Some Comments:

1) If  \alpha is negative the system is causal, and the filtered version g(t) of the signal f(t) will be causal.

2) Computing the smoothed signal does not save any computations with regards to computing the Laplace transform.

3) The derivation seems to show that their is a relationship between Laplace trancform and a windowed Fouier transform of the filtered signal.

4) To compute the Laplace transform based on the orginal signal use equation (5). To compute it based on the filtered signal use equation (11).

August 30, 2009 Posted by | Math | 8 Comments

Coriolis Forces

A derivation  for coriolis forces can be found be found on the wikipedia page for fictitious forces.

In general for an accelerating reference frame in rectangular coordinates the factious forces are given by:

\mathbf{F}_{\mbox{fictitious}} =-m\ \mathbf{a}_{AB} -2m\ \sum_{j=1}^3 v_j \ \frac{d \mathbf{u}_j}{dt} - m\ \sum_{j=1}^3 x_j \ \frac{d^2 \mathbf{u}_j}{dt^2}\ .


he first term is the Coriolis force, the second term is the centrifugal force, and the third term is the Euler force. When the rate of rotation doesn’t change, as is typically the case for a planet, the Euler force is zero.

Looking specifically at the Coriolis force:

- 2 m \boldsymbol \omega  \times \mathbf{v}

which gives in (east-west,north-sourth, height) coordinates:

\boldsymbol{ \Omega} = \omega \begin{pmatrix} 0 \\ \cos \varphi \\ \sin \varphi \end{pmatrix}\ , \boldsymbol{ v} = \begin{pmatrix} v_e \\ v_n \\ v_u \end{pmatrix}\ ,

\boldsymbol{ a}_C =-2\boldsymbol{\Omega \times v}= 2\,\omega\, \begin{pmatrix} v_n \sin \varphi-v_u \cos \varphi \\ -v_e \sin \varphi \\ v_e \cos\varphi\end{pmatrix}\ .

Where \phi is the latitudinal coordinate (equator=zero latitude).

In general the cross product for a coordinate with orthonormal direction vectors is given by:

\boldsymbol{\Omega \times v} = \begin{vmatrix} \boldsymbol{i}&\boldsymbol{j}&\boldsymbol{k} \\ \Omega_x & \Omega_y & \Omega_z \\ v_x & v_y & v_z \end{vmatrix}\ = \begin{pmatrix} \Omega_y v_z - \Omega_z v_y \\ \Omega_z v_x - \Omega_x v_z \\ \Omega_x v_y - \Omega_y v_x \end{pmatrix}\ ,

since the basis direction vectors are orthogonal in hopkins and simmons coordinates write write:

\boldsymbol{\Omega \times v} = \begin{vmatrix} \boldsymbol{e_{\sigma}}&\boldsymbol{e_{\mu}}&\boldsymbol{e_{\lambda}} \\ \Omega_{\sigma} & \Omega_{\mu} & \Omega_{\lambda} \\ W & V & U \end{vmatrix}\ = \begin{pmatrix} \Omega_{\mu} U - \Omega_{\lambda} V \\ \Omega_{\lambda} W - \Omega_{\sigma} U \\ \Omega_{\sigma} V - \Omega_{\mu} W \end{pmatrix}\ ,

(note with regards to weather the system is right handed we can choose the direction of the logitudanal cordinate \lambda to make it right handed.)

Just to recall from the post (Hoskins and Simmons (1974) Coordinate System):

\mu = sin( \theta ) where theta is the latitude.
\sigma = pressure/P_* Where P_*  is the surface pressure and \sigma is the vertical coordinate.
\lambda is the longitude.


U is the longitudinal component of the velocity
V is the latitudinal component of the velocity
W is the vertical component of the veolicty (not used in Hopkins and Simmons 1974)

Now the angular velocity of the earth in Hopkins and Simmons is given by:

\boldsymbol{ \Omega} = \omega \begin{pmatrix} \cos \varphi \\ \sin \varphi \\ 0 \end{pmatrix}\ = \omega \begin{pmatrix} \mp \sqrt{1-\mu^2} \\ \mu \\ 0 \end{pmatrix}\, \boldsymbol{ v} = \begin{pmatrix} W \\ V \\ U \end{pmatrix}\ ,

Where the sign of \pm is positive for the northern hemisphere and negative for the southern hemisphere.


\boldsymbol{\Omega \times v} =  \begin{pmatrix} \mu U \\  \pm \sqrt{1-\mu^2} U \\ \mp \sqrt{1-\mu^2} V - \mu  W \end{pmatrix} \

Some comments:

The result obtained is essentialy the same result that one would get if, they took the(east west,north south, altitude) coordinate system and replaced \phi with asin(\mu) .

The only differences are the order and sign of the components.  These are the only differences because both coordinate contain the same unit vectors.  In my example of a Hopkins and Simon’s like coordinate system I used a different order for the components then was used in my example for the (east-west, north south altitude) coordinate system. This will effect the sign in the cross product.

I wrote the z component of the angular velocity as \mp \sqrt{1-\mu^2} to emphasize that the positive direction for the z component in Simpons coordinate system  is downward. However, the actual angular rotation of the earth in simons coordinate system still have a postive \mu component depending on the which direction is defined as positive for the longitudinal  coordinate.

The order which we specify the coordinates determines the right handedness of the coordinate system.  Therefore, righthandedness is not inherently a geometric property because it depends on the order of the coordinates. For instance, in standard Cartesian coordinates e_z \times e_y = -e_x

In our case the first coordinat,e \sigma , was specified in the downward direction, our second coordinate,  \mu , points south, now using the right hand  rule means that gives the positive direction for the third coordinate \lambda in the east direction.

It is for these reasons that differences can arrise, and therefore it is very important when doing cross products to clearly express the postive direction of the coordinate unit vectors and the order of the coordinates.

August 29, 2009 Posted by | GCM (General Circulation Models) | 2 Comments

Vector Operations in Hoskins and Simmons Coordinates

In my post Hoskins and Simmons (1974) Coordinate System, I derived the following scaling quantities which will be used to derive the vector operations of Grad, Div and Curl in Hoskins Coordinate system.

h_{\mu}= \left| \frac{r}{\sqrt{1-\mu^2}} \right|
h_{\lambda}=r \mu
h_{\sigma} =-\sigma \frac{mg}{RT}

The coordinates in Hoskins coordinate system are dimensionless . (see nondimentionalization of Navier Stokes).

The gradient is defined as (see lectures on coordinate transforms):

\nabla \phi = \frac{\vec e_{\mu}}{h_{\mu}} \frac{\partial \phi}{\partial \mu} + \frac{\vec e_{\lambda}}{h_{\lambda}}\frac{\partial \phi}{\partial \lambda}+\frac{\vec e_{\sigma}}{h_{\sigma}}\frac{\partial \phi}{\partial \sigma}=\left(\vec e_{\mu}\right)\frac{1-\mu^2}{r} \frac{\partial \phi}{\partial \mu} + \left(\vec e_{\lambda}\right)\frac{1}{r \mu}\frac{\partial \phi}{\partial \lambda}-\left(\vec e_{\sigma}\right)\frac{RT}{\sigma mg}\frac{\partial \phi}{\partial \sigma}

The divergence is defined as:

\nabla \cdot = \frac{1}{h_{\mu}} \frac{\partial }{\partial \mu} + \frac{1}{h_{\lambda}}\frac{\partial}{\partial \lambda}+\frac{1}{h_{\sigma}}\frac{\partial}{\partial \sigma}=\frac{1-\mu^2}{r} \frac{\partial}{\partial \mu} + \frac{1}{r \mu}\frac{\partial}{\partial \lambda}-\frac{RT}{\sigma mg}\frac{\partial}{\partial \sigma}

The curl is defined by:

\nabla \times \vec A = \frac{1}{h_{\sigma}h_{\mu}h_{\lambda}} \left| \begin{array}{ccc} \vec e_{\sigma} & \vec e_{\mu} & \vec e_{\lambda} \\ \frac{\partial}{\partial \sigma} & \frac{\partial}{\partial \mu} & \frac{\partial}{\partial \lambda} \\ A_{\sigma}h_{\sigma} & A_{\mu}h_{\mu} &  A_{\lambda}h_{\lambda} \end{array} \right|

(note, the direction of the longitudinal coordinate \lambda is defined to obey the right hand rule)

This gives for the components

\left[ \nabla \times \vec A \right]_{\sigma}=\frac{RT}{r \mu \sigma mg}\left(\sigma \frac{mg}{RT}\frac{\partial A_{\lambda}}{\partial \mu}+ r \mu\frac{\partial A_{\mu}}{\partial \lambda} \right)

\left[ \nabla \times \vec A \right]_{\mu}= -\frac{\sqrt{1-\mu^2}RT}{r\sigma mg}\left( \frac{r}{\sqrt{1-\mu^2}}\frac{\partial A_{\sigma}}{\partial \lambda}+\sigma \frac{mg}{RT} \frac{\partial A_{\lambda}}{\partial \sigma} \right)

\left[ \nabla \times \vec A \right]_{\lambda}= \frac{\sqrt{1-\mu^2}RT}{r^2 \mu }\left( r \mu\frac{\partial A_{\mu}}{\partial \sigma}-\frac{r}{\sqrt{1-\mu^2}} \frac{\partial A_{\sigma}}{\partial \mu} \right)

Which Simplifies to:

\left[ \nabla \times \vec A \right]_{\sigma}= \frac{1}{r\mu}\frac{\partial A_{\lambda}}{\partial \mu}+ \frac{RT}{\sigma m g}\frac{\partial A_{\mu}}{\partial \lambda}

\left[ \nabla \times \vec A \right]_{\mu}=  -\frac{RT}{\sigma m g}\frac{\partial A_{\sigma}}{\partial \lambda}- \frac{\sqrt{1-\mu^2}}{r} \frac{\partial A_{\lambda}}{\partial \sigma}

\left[ \nabla \times \vec A \right]_{\lambda}=  \frac{\sqrt{1-\mu^2}RT}{r}\frac{\partial A_{\mu}}{\partial \sigma}-\frac{RT}{r\mu} \frac{\partial A_{\sigma}}{\partial \mu}

August 29, 2009 Posted by | Uncategorized | 1 Comment

API/Object Viewers/Memory Mapping/

The more code a programmer can reuse the more efficient they can be. In windows this could mean reusing com/ole components and other APIs. Here are two useful programs for viewing APIs:

OLE/COM Object Explorer 1.1

Windows API Viewer

I was inquiring about how to manage the transfer of large amounts of data between programs and I was pointed to two interesting concepts:

All modern operating systems include a facility called “memory mapping,” which maps a range of addresses in the program’s virtual address space to a file. If you read from those addresses, you’ll get data from the file. It is up to the operating system to determine whether to load the data into RAM all at once, or to read it from the disk in chunks as necessary.
If you’re trying to share large amounts of memory between two programs running on the same computer, you should note that all modern operating systems provide mechanisms for shared memory. These shared memory segments can be mapped into the virtual address space of multiple programs simultaneously. Two or more programs can read or write to the shared memory exactly as if it were normal, private memory. (But you should include some thread-safety mechanisms, like mutexes, to make sure your programs won’t step on each other’s toes.)

If you’re trying to share large amounts of memory between programs running on separate computers, use MPI or some other multi-processing library.

Here is what wikipedia has to say about memory maps:

The primary benefit of memory mapping a file is increased I/O performance, especially when used on small files. Accessing memory mapped files is faster than using direct read and write operations for two reasons. Firstly, a system call is orders of magnitude slower than a simple change of program’s local memory. Secondly, in most operating systems the memory region mapped actually is the kernel’s file cache, meaning that no copies need to be created in user space. Using system calls would inevitably involve the time consuming operation of memory copying.

Certain application level memory-mapped file operations also perform better than their physical file counterparts. Applications can access and update data in the file directly and in-place, as opposed to seeking from the start of the file or rewriting the entire edited contents to a temporary location. Since the memory-mapped file is handled internally in pages, linear file access (as seen, for example, in flat file data storage or configuration files) requires disk access only when a new page boundary is crossed, and can write larger sections of the file to disk in a single operation.

A possible benefit of memory-mapped files is a “lazy loading”, thus using small amounts of RAM even for a very large file. Trying to load the entire contents of a file that is significantly larger than the amount of memory available can cause severe thrashing as the operating system reads from disk into memory and simultaneously pages from memory back to disk. Memory-mapping may not only bypass the page file completely, but the system only needs to load the smaller page-sized sections as data is being edited, similarly to demand paging scheme used for programs.

The windows utility to do this is called CreateFileMapping

As for multi-processing library the following was recommended:

I haven’t found much but the following wikipedia link seems relevant:

August 29, 2009 Posted by | Computer Science and Modeling | 2 Comments