It looks like you're new here. If you want to get involved, click one of these buttons!

- All Categories 2.2K
- Applied Category Theory Course 356
- Applied Category Theory Seminar 4
- Exercises 149
- Discussion Groups 50
- How to Use MathJax 15
- Chat 482
- Azimuth Code Project 108
- News and Information 145
- Azimuth Blog 149
- Azimuth Forum 29
- Azimuth Project 189
- - Strategy 108
- - Conventions and Policies 21
- - Questions 43
- Azimuth Wiki 711
- - Latest Changes 701
- - - Action 14
- - - Biodiversity 8
- - - Books 2
- - - Carbon 9
- - - Computational methods 38
- - - Climate 53
- - - Earth science 23
- - - Ecology 43
- - - Energy 29
- - - Experiments 30
- - - Geoengineering 0
- - - Mathematical methods 69
- - - Meta 9
- - - Methodology 16
- - - Natural resources 7
- - - Oceans 4
- - - Organizations 34
- - - People 6
- - - Publishing 4
- - - Reports 3
- - - Software 21
- - - Statistical methods 2
- - - Sustainability 4
- - - Things to do 2
- - - Visualisation 1
- General 41

Options

I will start with this tomorrow after lunch (here) John. And I decided to have a separate chat thread for that! I'll read the paper first, so you will probably see questions after that.

It's probably good to start with the simplest models and work our way up. Earlier here on the Forum we had a lot of fun discussing models that exhibit abrupt climate change in response to changes in the "insolation" (very roughly, the amount of sunlight hitting the Earth). Some of these models have been proposed to study the "Snowball Earth" era about 850 million years ago, and also the more recent glacial cycles of the last 1.8 million years or so.

In order of increasing complexity, these are four simple models:

1) In our page Snowball Earth we discuss a paper by I. Zaliapin and M. Ghil called "Another look at climate sensitivity". There's a link to this paper. Equations (11), (12), (13) and (14) in this paper describe a simple energy balance model. There's really just a single equation in the end, equation (15). If you unravel all the notation, it says:

$$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ where $T$ is temperature as a function of time $t$, and $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$ are constants.

## Comments

Great!

`Great!`

How do you want to do 1) John? My intention was to use Sage and get the relevant equations and start to see what I can do both with analytical and numerical features of Sage. Some Qs

I assume you want to focus on slow dynamics in 3. which they also do in the paper?

Are you thinking of using the same approach as you did with the Hopf bifurcation, but with Saddle-node instead?

If so how would you like to add noise?

`How do you want to do 1) John? My intention was to use Sage and get the relevant equations and start to see what I can do both with analytical and numerical features of Sage. Some Qs I assume you want to focus on slow dynamics in 3. which they also do in the paper? Are you thinking of using the same approach as you did with the Hopf bifurcation, but with Saddle-node instead? If so how would you like to add noise?`

For me the interesting story is that there is (or could be) more than one meta-stable state and that there is another locally stable state that is radically different from the one observed today. That is interesting because a lot of the discussion of climate change revolves around the statements that

a) climate has always been changing,

b) but it has been stable within reasonable bounds for millions of years.

The model can show that the combination of a trend and of stochastic forcing can easily push the climate over the barrier from one locally stable state to another, but does so using a physically reasonable, if very oversimplified, approximation to the thermodynamics of the atmosphere, unlike the generic example of stochastic resonance.

`For me the interesting story is that there is (or could be) more than one meta-stable state and that there is another locally stable state that is radically different from the one observed today. That is interesting because a lot of the discussion of climate change revolves around the statements that a) climate has always been changing, b) but it has been stable within reasonable bounds for millions of years. The model can show that the combination of a trend and of stochastic forcing can easily push the climate over the barrier from one locally stable state to another, but does so using a physically reasonable, if very oversimplified, approximation to the thermodynamics of the atmosphere, unlike the generic example of [[stochastic resonance]].`

Yes I assume that it is feasible to use $\mu$ and follow how it is done in the paper and then gradually add noise in some form that is also motivated as a physical property. So it is just to find the correct ones

`Yes I assume that it is feasible to use $\mu$ and follow how it is done in the paper and then gradually add noise in some form that is also motivated as a physical property. So it is just to find the correct ones`

Staffan wrote:

What do you you mean by "3."?

In the Zaliapin-Ghil paper they take this equation

$$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ where $T$ is temperature as a function of time $t$, and $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$ are constants. Then in section 3.3.2 they consider the "slow, quasi-adiabatic" dynamics of this equation. That's just a fancy of way of saying they assume

$$ \frac{d T}{d t} = 0$$ and solve for $T$ as a function of $\mu$, holding other things fixed at specific values. Here $\mu Q_0$ is the '"insolation", that is, very roughly, the amount of sunlight reaching the Earth. $Q_0$ is some "standard" amount of insolation, and $\mu$ is some parameter near $1$ which they call the "fractional insolation change". They get this graph:

So, we see that if $\mu$ is low, there is just one solution for $T$, which we could call the "Snowball Earth" solution. As we increase $\mu$ we reach a value there there are three solutions: two stable ones drawn as solid curves (the "Snowball Earth" solution and the "Hot Earth" solution) together with an intermediate unstable one drawn as a dashed curve. As we increase $\mu$ even more there is only one solution: the "Hot Earth" solution.

So far this is not exciting enough to deserve a computer program, but let me stop here for now.

`Staffan wrote: > I assume you want to focus on slow dynamics in 3. which they also do in the paper? What do you you mean by "3."? In the [Zaliapin-Ghil paper](http://www.nonlin-processes-geophys.net/17/113/2010/npg-17-113-2010.pdf) they take this equation $$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ where $T$ is temperature as a function of time $t$, and $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$ are constants. Then in section 3.3.2 they consider the "slow, quasi-adiabatic" dynamics of this equation. That's just a fancy of way of saying they assume $$ \frac{d T}{d t} = 0$$ and solve for $T$ as a function of $\mu$, holding other things fixed at specific values. Here $\mu Q_0$ is the '"insolation", that is, very roughly, the amount of sunlight reaching the Earth. $Q_0$ is some "standard" amount of insolation, and $\mu$ is some parameter near $1$ which they call the "fractional insolation change". They get this graph: <img src = "http://math.ucr.edu/home/baez/zaliapin_ghil_snowball_earth.jpg" alt = ""/> So, we see that if $\mu$ is low, there is just one solution for $T$, which we could call the "Snowball Earth" solution. As we increase $\mu$ we reach a value there there are three solutions: two stable ones drawn as solid curves (the "Snowball Earth" solution and the "Hot Earth" solution) together with an intermediate unstable one drawn as a dashed curve. As we increase $\mu$ even more there is only one solution: the "Hot Earth" solution. So far this is not exciting enough to deserve a computer program, but let me stop here for now.`

Next, we might take the same equation:

$$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ but treat $\mu$ not as a constant but as varying as a function of time. The simplest interesting case would be to let $\mu$ oscillate around an average value $\mu_0$:

$$\mu(t) = \mu_0 + a \cos(\omega t) $$ It'll be interesting if the oscillations are big enough that:

when $\mu$ is large the only stable solution (from the "adiabatic" analysis in my last comment) is the "Hot Earth" solution, but

when $\mu$ is small the only stable solution is the "Snowball Earth" solution.

Then if we start with $T(0)$ "hot" and the oscillation frequency $\omega$ is not too fast, the temperature $T(t)$ should stay "Hot" until $\mu$ gets sufficiently small, at which point the Earth should "snap" rather suddenly to the "Snowball" solution... and then when $\mu$ gets big enough it should "snap back" to the "Hot" solution.

The next step would be to add noise, but it would be nice to see this first. In particular, it would be nice if we could agree on some interesting values of all these parameters $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$, and see what this graph looks like for those parameters:

That'll tell us interesting values for the parameters $\mu_0$ and $a$ in this formula:

$$\mu(t) = \mu_0 + a \cos(\omega t) $$ Of course all of this is just an oversimplified model - not much to do with the actual Earth, except in a rough way. Things get a bit more realistic when we reach the Crucifix-Rougier model.

`Next, we might take the same equation: $$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ but treat $\mu$ not as a constant but as varying as a function of time. The simplest interesting case would be to let $\mu$ oscillate around an average value $\mu_0$: $$\mu(t) = \mu_0 + a \cos(\omega t) $$ It'll be interesting if the oscillations are big enough that: * when $\mu$ is large the only stable solution (from the "adiabatic" analysis in my last comment) is the "Hot Earth" solution, but * when $\mu$ is small the only stable solution is the "Snowball Earth" solution. Then if we start with $T(0)$ "hot" and the oscillation frequency $\omega$ is not too fast, the temperature $T(t)$ should stay "Hot" until $\mu$ gets sufficiently small, at which point the Earth should "snap" rather suddenly to the "Snowball" solution... and then when $\mu$ gets big enough it should "snap back" to the "Hot" solution. The next step would be to add noise, but it would be nice to see this first. In particular, it would be nice if we could agree on some interesting values of all these parameters $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$, and see what this graph looks like for those parameters: <img src = "http://math.ucr.edu/home/baez/zaliapin_ghil_snowball_earth.jpg" alt = ""/> That'll tell us interesting values for the parameters $\mu_0$ and $a$ in this formula: $$\mu(t) = \mu_0 + a \cos(\omega t) $$ Of course all of this is just an oversimplified model - not much to do with the actual Earth, except in a rough way. Things get a bit more realistic when we reach the [Crucifix-Rougier model](http://www.azimuthproject.org/azimuth/show/Bayesian%20prediction%20of%20the%20next%20glacial%20inception).`

Sorry I meant section 3, like you talked about here. I got the basic equation into Sage by just copying the source version here and pasting into a Sage cell and adding a keyword %latex to the cell (This is just for me to remember. I'll continue a bit and add time varying $\mu$ . Should we have page here on Azimuth now or should we wait?

`Sorry I meant section 3, like you talked about here. I got the basic equation into Sage by just copying the source version here and pasting into a Sage cell and adding a keyword %latex to the cell (This is just for me to remember. I'll continue a bit and add time varying $\mu$ . Should we have page here on Azimuth now or should we wait?`

Staffan wrote:

Okay, I was afraid you might have meant stage 3) in this grand plan. I'm glad you're still working on 1) and 2) here:

`Staffan wrote: > Sorry I meant section 3, like you talked about here. Okay, I was afraid you might have meant stage 3) in this grand plan. I'm glad you're still working on 1) and 2) here: > In order of increasing complexity, these are four simple models: > 1) In our page [[Snowball Earth]] we discuss a paper by I. Zaliapin and M. Ghil called "Another look at climate sensitivity". There's a link to this paper. Equations (11), (12), (13) and (14) in this paper describe a simple energy balance model. There's really just a single equation in the end, equation (15). If you unravel all the notation, it says: > $$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh \kappa(T - T_c)}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ > where $T$ is temperature as a function of time $t$, and $c, Q_0, c_1, c_2, \kappa, T_c, \sigma, m, T_0$ are constants. > 2) I can also imagine a stochastic version of the Zaliapin-Ghil model, which should be more interesting, but more complicated. > 3) A more complicated deterministic model is described in [[Bayesian prediction of the next glacial inception]]: > $$ \frac{d}{d t}I(t) = -a_1 \left(k_{\mu} \mu(t) + k_{\theta} \theta(t) - K_I I(t)\right) + k_R R(t) $$ > $$ \frac{d}{d t}\mu(t) = b_1 \mu(t) - b_2 \mu(t)^2 - b_3 \mu(t)^3 - b_{\theta} \theta(t) $$ > $$ \frac{d}{d t}\theta(t) = - c_1 I(t) - K_{\theta} \theta(t) $$ > The meaning of these equation is described on that page, and in more detail in the paper by Crucifix and Rougier, which the page links to. That paper also proposes values for all the constants involved here. > 4) The stochastic version of the Crucifix-Rougier model is even more interesting and complicated. That's what their paper is really about.`

There is a small typo in your unraveled eq.17 the ice albedo feedback should be:

$$c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2})$$

`There is a small typo in your unraveled eq.17 the ice albedo feedback should be: $$c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2})$$`

I left the Sage draft worksheet public, so you can see and comment: I have some small things to "unravel" in the code but you might start thinking about what to try for parameter values, or double checking constants in the Sage code with zhalian,ghil paper. But you will not see the @interactive part unless you register at:

this alternative server. It is not so heavily loaded, so I have been able to work more concentrated :-)

`I left the Sage draft worksheet public, so [you can see](http://demo.sagenb.org/home/pub/92/) and comment: I have some small things to "unravel" in the code but you might start thinking about what to try for parameter values, or double checking constants in the Sage code with zhalian,ghil paper. But you will not see the @interactive part unless you register at: [this alternative server](http://demo.sagenb.org). It is not so heavily loaded, so I have been able to work more concentrated :-)`

Has anyone verified Figure 5? I tried, but failed.

It's roughly the right shape, but looks off by a factor of 2 (or worse). Here is the code:

I'm new to Python so could be doing something silly.

An idea I have to take this in the desired direction of a dynamic insolation driver is to

PS: Only a mathematician would be so unthoughtful about notation as to write down something like $\mu(0) = \mu_0 + a$ ;)

`Has anyone verified Figure 5? I tried, but failed. <img src="http://www.azimuthproject.org/azimuth/files/ZG.png" width = "500" alt = ""/> It's roughly the right shape, but looks off by a factor of 2 (or worse). Here is the code: import numpy import scipy import matplotlib from pylab import * Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0 = (1.9e-15)**(-6) sigma = 5.6697e-8; Tc = 273.0; kappa= 0.05 T = arange(100.0, 400.0+1.0, 1.0) mu = (sigma*T**4*(1 - m*tanh((T/T0)**6)))/(Q0*(1 - (c1 + .5*c2*(1 - tanh(kappa*(T - Tc)))))) plot(mu, T, '-', lw=2) xlabel('Fractional Insolation Change (mu)') ylabel('Temperature (T)') grid(True) show() I'm new to Python so could be doing something silly. An idea I have to take this in the desired direction of a dynamic insolation driver is to 1. Get the steady state correct (which means getting this chart right). 1. Begin at some point on the chart that is in steady state. This will be our initial conditions. 1. Turn on a driver $\mu(t) = \mu_0 + u(t) a sin(\omega t)$, where $u(t)$ is the unit step function. 1. Solve the time evolution of $T$ numerically. PS: Only a mathematician would be so unthoughtful about notation as to write down something like $\mu(0) = \mu_0 + a$ ;)`

I'm new here; I hope I'll be able to find the time to contribute more than just typo-finding.

`Eric: your factor of 2 comes from a small typo: "T0 = (1.9e-15)**(-6)" should be "T0 = (1.9e-15)**(1/6.)". I'm new here; I hope I'll be able to find the time to contribute more than just typo-finding.`

Thanks Matt! That's probably it. Or maybe **(-1/6). I'll fix it when I get home and report back.

`Thanks Matt! That's probably it. Or maybe **(-1/6). I'll fix it when I get home and report back.`

Yep. That was definitely a typo, but it still doesn't fix things completely. The new code looks like:

One interesting difference between Matlab and Python is that you need to pay attention to data types. For example, in Python

comes out to be

It took a second to figure out why. It is because Python seems to be rounding -1/6 down to the nearest integer, which is -1. So it is computing (1.9e-15)**(-1).

Here's the new chart:

Still not quite right. I'm missing something obvious...

`Yep. That was definitely a typo, but it still doesn't fix things completely. The new code looks like: import numpy import scipy import matplotlib from pylab import * Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0 = (1.9e-15)**(-1.0/6.0) sigma = 5.6697e-8; Tc = 273.0; kappa= 0.05 T = arange(100.0, 400.0+1.0, 1.0) mu = (sigma*T**4*(1 - m*tanh((T/T0)**6)))/(Q0*(1 - (c1 + .5*c2*(1 - tanh(kappa*(T - Tc)))))) plot(mu, T, '-', lw=2) xlabel('Fractional Insolation Change (mu)') ylabel('Temperature (T)') grid(True) show() One interesting difference between Matlab and Python is that you need to pay attention to data types. For example, in Python T0 = (1.9e-15)**(-1/6) comes out to be T0 = 526315789473684.2 It took a second to figure out why. It is because Python seems to be rounding -1/6 down to the nearest integer, which is -1. So it is computing (1.9e-15)**(-1). Here's the new chart: <img src="http://www.azimuthproject.org/azimuth/files/ZG_II.png" width = "500" alt = ""/> Still not quite right. I'm missing something obvious...`

You're hitting a Python version specific thing. Before version 3.0, Python decided to copy the C behaviour that int/int= result rounded down to nearest below int. For version 3 onwards it was decided "the C behaviour is counter-to-reasonable-expectations, so we'll change it to int/int=float". (Incidentally, it appears Java has the C behaviour, so beware of this when doing anything in Java as well.)

(I can't comment on any issues in the rest of the code because I haven't been following the original paper you're reproducing.)

`You're hitting a Python version specific thing. Before version 3.0, Python decided to copy the C behaviour that int/int= result rounded down to nearest below int. For version 3 onwards it was decided "the C behaviour is counter-to-reasonable-expectations, so we'll change it to int/int=float". (Incidentally, it appears Java has the C behaviour, so beware of this when doing anything in Java as well.) (I can't comment on any issues in the rest of the code because I haven't been following the original paper you're reproducing.)`

Hmm. My own sign mistake (i.e. 1/6 in the exponent rather than -1/6) actually reproduces their figure. So I think they messed up: they flipped the sign, giving essentially a flat 0.6 for $1 - m \tanh((T/T_0)^6)$ rather than a function that interpolates between $1$ at $T \to 0$ and $0.6$ for $T \gt 350 \mathrm{K}$. Nothing qualitative changes, but I think quantitatively the plot in their paper is wrong and yours is now correct.

`Hmm. My own sign mistake (i.e. 1/6 in the exponent rather than -1/6) actually reproduces their figure. So I think they messed up: they flipped the sign, giving essentially a flat 0.6 for $1 - m \tanh((T/T_0)^6)$ rather than a function that interpolates between $1$ at $T \to 0$ and $0.6$ for $T \gt 350 \mathrm{K}$. Nothing qualitative changes, but I think quantitatively the plot in their paper is wrong and yours is now correct.`

Eric,

I reproduce your figure when I implement the same equations in R. I tried to implement it right from the paper, without copying your code (except for the constants, which I double checked, including going back to the literature to verify the T0 value).

`Eric, I reproduce your figure when I implement the same equations in R. I tried to implement it right from the paper, without copying your code (except for the constants, which I double checked, including going back to the literature to verify the T0 value). Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0 = (1.9e-15)^(-1/6) sigma = 5.6697e-8; Tc = 273; kappa = 0.05 T = seq(100, 400, by=1) alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2 lhs = Q0*(1-alpha) g = 1 - m*tanh((T/T0)^6) rhs = sigma*T^4 * g mu = rhs/lhs # mu Q0(1-alpha) = sigma g T^4 plot(mu, T, xlim=c(0,3), ylim=c(100,400), type="l", lwd=2, col="blue", xlab="Fractional insolation change (mu)", ylab="Temperature (T)", main="Zaliapin and Ghil (2010), Fig. 5", xaxs="i", yaxs="i") abline(h=seq(100, 400, by=50), lty="dotted", col="gray") abline(v=seq(0, 2.5, by=0.5), lty="dotted", col="gray")`

Matt,

I think that's it. When plotting they used 1/6 instead of -1/6 defining T0. Their equations themselves are correct.

`Matt, I think that's it. When plotting they used 1/6 instead of -1/6 defining T0. Their equations themselves are correct.`

Updated R code:

`Updated R code: Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4 sigma = 5.6697e-8; Tc = 273; kappa = 0.05 T0 = (1.9e-15)^(-1/6) # correct ##T0 = (1.9e-15)^(1/6) # ZG10 bug T = seq(100, 400, by=1) # fractional insolation change mu = function(T) { alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2 lhs = Q0*(1-alpha) g = 1 - m*tanh((T/T0)^6) rhs = sigma*T^4 * g mu = rhs/lhs # mu Q0(1-alpha) = sigma g T^4 } # EBM bifurcation diagram plot(mu(T), T, xlim=c(0,3), ylim=c(100,400), type="l", lwd=2, col="blue", xlab="Fractional insolation change (mu)", ylab="Temperature (T)", main="Zaliapin and Ghil (2010), Fig. 5 (corrected)", xaxs="i", yaxs="i") abline(h=seq(100, 400, by=50), lty="dotted", col="gray") abline(v=seq(0, 2.5, by=0.5), lty="dotted", col="gray") # modern Earth state T = 300 points(mu(T), T, pch=21, col="black", bg="green", cex=2)`

Ah. Thanks for confirming :)

Now the fun can begin :)

`Ah. Thanks for confirming :) Now the fun can begin :)`

I just used their bug value $T_0^{-6}$ but took care to not raise it ZG eq (14) and got the correct plots for figure ZG fig 4, so I went straight for next step to add time varying mu, see link above to open nb. But I will double check and plot fig 5 in Sage as well . I added Nathan's R code to the Sage notebook. But if You want to run that R cell in the notebook you need to either run it on your own machine with Sage installed or change the R code to save instead of plotting.

`I just used their bug value $T_0^{-6}$ but took care to not raise it ZG eq (14) and got the correct plots for figure ZG fig 4, so I went straight for next step to add time varying mu, see link above to open nb. But I will double check and plot fig 5 in Sage as well . I added Nathan's R code to [the Sage notebook](http://demo.sagenb.org/home/pub/92/). But if You want to run that R cell in the notebook you need to either run it on your own machine with Sage installed or change the R code to save instead of plotting.`

# steady state in zg and their constant values from 3.2

Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0_6= 1.9e-15

sigma_sb= 5.6697e-8; Tc = 273.0

plot zg figure 4 to see if correct

var('T')

mu = 1.0; kappa= 0.05

zg16 = muQ0(1 - (c1 + c2(1 - tanh(kappa(T - Tc)))/2.0)) == sigma_sbT^4 *(1 - m tanh((T^6/T0_6)))

zg16ri = zg16.lhs()

zg16ro = zg16.rhs()

pro = plot(zg16ro,(T,150,400),rgbcolor='red')

pri = plot(zg16ri,(T,150,400))

add the two plots to see the 3 fixed points when mu= 1

pri + pro

`Here is that sage code - note it knows about symbolic lhs, rhs of an equation, zg16 and it speaks that universal language that starts with an m :-) # steady state in zg and their constant values from 3.2 Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4; T0_6= 1.9e-15 sigma_sb= 5.6697e-8; Tc = 273.0 # plot zg figure 4 to see if correct var('T') mu = 1.0; kappa= 0.05 zg16 = mu*Q0*(1 - (c1 + c2*(1 - tanh(kappa*(T - Tc)))/2.0)) == sigma_sb*T^4 *(1 - m* tanh((T^6/T0_6))) zg16ri = zg16.lhs() zg16ro = zg16.rhs() pro = plot(zg16ro,(T,150,400),rgbcolor='red') pri = plot(zg16ri,(T,150,400)) # add the two plots to see the 3 fixed points when mu= 1 pri + pro`

My code shouldn't have 'border=TRUE' in the last line; I forgot to take that out. I edited the code above accordingly.

`My code shouldn't have 'border=TRUE' in the last line; I forgot to take that out. I edited the code above accordingly.`

Ok I'll changed that in that in the Sage worksheet. I also changed my Sage equations and constants so they agree with everybody else here who uses $tanh((\frac {T} {T_0})^6)$

`Ok I'll changed that in that in the Sage worksheet. I also changed my Sage equations and constants so they agree with everybody else here who uses $tanh((\frac {T} {T_0})^6)$`

...and when i try this i get the plots in fig4 wrong . when kappa= 0.05 and mu = 1 and T_0 with your value you should get three fix points, but I only get one. Nathan did you also plot fig4 ?

`...and when i try this i get the plots in fig4 wrong . when kappa= 0.05 and mu = 1 and T_0 with your value you should get three fix points, but I only get one. Nathan did you also plot fig4 ?`

They have the same $T_0$ bug in Figure 4. With that bug, I replicate their figure.

`They have the same $T_0$ bug in Figure 4. With that bug, I replicate their figure. Q0 = 342.5; c1 = 0.15; c2 = 0.7; m = 0.4 sigma = 5.6697e-8; Tc = 273; kappa = 0.05 T0 = (1.9e-15)^(-1/6) # correct ##T0 = (1.9e-15)^(1/6) # ZG10 bug T = seq(150, 400, by=1) # incoming radiation Ri = function(T, mu) { alpha = c1 + c2*(1-tanh(kappa*(T-Tc)))/2 Ri = mu * Q0 * (1-alpha) } # outgoing radiation Ro = function(T) { g = 1 - m*tanh((T/T0)^6) Ro = sigma * g * T^4 } # plot outgoing radation as a function of temperature plot(T, Ro(T), xlim=c(150,400), ylim=c(0,630), type="l", lty="dashed", lwd=2, col="darkgreen", xlab="Temperature, T", ylab="Radiation, R", main="Zaliapin and Ghil (2010), Fig. 4", xaxs="i", yaxs="i") # plot incoming radiation, for different fractional insolations lines(T, Ri(T, mu=0.5), lwd=2, col="blue") lines(T, Ri(T, mu=1.0), lwd=2, col="blue") lines(T, Ri(T, mu=2.0), lwd=2, col="blue")`

ok, thanks that explains things . I tried both your R scripts in my local R installaton and you are correct. Its a "science bug"!

I also added the other R script to the Sage worksheet.

`ok, thanks that explains things . I tried both your R scripts in my local R installaton and you are correct. Its a "science bug"! I also added the other R script to the Sage worksheet.`

Eric: Have you found any fun parameters for the added mu oscillation ? In the worksheet I intend to look at a two things: just having a fixed sine and adjusting that with one parameter - strength which applies to the whole $\mu$ in ZG10 to see what happens (first I was doing an interactive version which is still there, but I have a feeling that it becomes to tedious to look for interesting things).

The other was too look into more realistic positive feedback oscillations that could potentially trigger a fold bifurcation into instability eg global homoclinic bifurcations,phase locking, where some examples would be anything simple that is incrementally changing in nature and all the way to natural disasters.

So should we do both or just settle for the first for now?

`Eric: Have you found any fun parameters for the added mu oscillation ? In the worksheet I intend to look at a two things: just having a fixed sine and adjusting that with one parameter - strength which applies to the whole $\mu$ in ZG10 to see what happens (first I was doing an interactive version which is still there, but I have a feeling that it becomes to tedious to look for interesting things). The other was too look into more realistic positive feedback oscillations that could potentially trigger a fold bifurcation into instability eg global homoclinic bifurcations,phase locking, where some examples would be anything simple that is incrementally changing in nature and all the way to natural disasters. So should we do both or just settle for the first for now?`

Hi Staffan,

I haven't had much free time lately, so haven't made any progress beyond what is here. I am a little bummed about their error though. It does seem to be a bit of a blow to a significant portion of their discussion.

Once I get some time, I planned to implement a finite difference (with central differencing) solution for a sinusoidal insolation just to play around. I would not be bummed if you beat me to it.

`Hi Staffan, I haven't had much free time lately, so haven't made any progress beyond what is here. I am a little bummed about their error though. It does seem to be a bit of a blow to a significant portion of their discussion. Once I get some time, I planned to implement a finite difference (with central differencing) solution for a sinusoidal insolation just to play around. I would not be bummed if you beat me to it.`

Don't worry I am also bummed about their error. But right now I left the Sage worksheet, with some interactive 3d plots so for now I leave it like that and ill fiddle around a bit now and then until John comes and spank my hands :-) I have some other ideas i want to put up on Azimuth before that.

`Don't worry I am also bummed about their error. But right now I left the Sage worksheet, with some interactive 3d plots so for now I leave it like that and ill fiddle around a bit now and then until John comes and spank my hands :-) I have some other ideas i want to put up on Azimuth before that.`

Hi guys! I'm sorry that I've been neglecting this climate project... I've been absorbed in that "network theory" project. They're both part of Azimuth and they're even both part of the same huge mega-project of understanding models of stochastic systems... but I should pay attention to you folks.

Can somebody summarize the progress so far? I'll try, and you can correct me:

1) You guys took the equation

$$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ and looked for equilibrium states with $d T / d t = 0$, getting

$$\mu = \frac{ \sigma T^4 (1 - m tanh((T/T_0)^6))}{Q_0 (1 - (c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2}))}$$ which expresses $\mu$ (fractional insolation change) as a function of $T$ (the Earth's average temperature).

2) You graphed $\mu$ as a function of $T$ and discovered that Zaliapin and Ghil's graph was wrong.

3) You even figured out what they did wrong: they used $T_0 = (1.9 \cdot 10^{-15})^{1/6}$ when creating their graph, instead of what they said they were using, $T_0 = (1.9 \cdot 10^{-15})^{-1/6}$.

4) Eric and Staffan are "bummed" because they're wrong.

All this sounds great except 4). You should be

happybecause they're wrong - and you caught it! Where's your scientific spirit? There's nothing more fun than finding a mistake in someone else's work! Iloveit when that happens!Besides, this will be a great advertisement for the Azimuth Project, and the whole idea of "open source science". Here we are, a bunch of "amateurs", and we start trying to replicate the simplest model we can find... and then, right away, you guys discover a mistake in the published literature!

Unfortunately it sounds like the mistake is too small to matter much... but still, it's amusing.

What next? Is there anything I can do to help?

I could help polish the page Another look at climate sensitivity, but I don't understand the graphs at the bottom: which graph came from which procedure. If someone writes a quick note to make that clear, I can add more detail.

`Hi guys! I'm sorry that I've been neglecting this climate project... I've been absorbed in that "network theory" project. They're both part of Azimuth and they're even both part of the same huge mega-project of understanding models of stochastic systems... but I should pay attention to you folks. Can somebody summarize the progress so far? I'll try, and you can correct me: 1) You guys took the equation $$ c \frac{d T}{d t} = \mu Q_0 (1 - (c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2})) - \sigma T^4 (1 - m tanh((T/T_0)^6)) $$ and looked for equilibrium states with $d T / d t = 0$, getting $$\mu = \frac{ \sigma T^4 (1 - m tanh((T/T_0)^6))}{Q_0 (1 - (c_1 + c_2 \frac{1 - tanh (\kappa(T - T_c))}{2}))}$$ which expresses $\mu$ (fractional insolation change) as a function of $T$ (the Earth's average temperature). 2) You graphed $\mu$ as a function of $T$ and discovered that [Zaliapin and Ghil's graph](http://www.azimuthproject.org/azimuth/show/Another+look+at+climate+sensitivity) was wrong. 3) You even figured out what they did wrong: they used $T_0 = (1.9 \cdot 10^{-15})^{1/6}$ when creating their graph, instead of what they said they were using, $T_0 = (1.9 \cdot 10^{-15})^{-1/6}$. 4) Eric and Staffan are "bummed" because they're wrong. All this sounds great except 4). You should be _happy_ because they're wrong - and you caught it! Where's your scientific spirit? There's nothing more fun than finding a mistake in someone else's work! I _love_ it when that happens! Besides, this will be a great advertisement for the Azimuth Project, and the whole idea of "open source science". Here we are, a bunch of "amateurs", and we start trying to replicate the simplest model we can find... and then, right away, you guys discover a mistake in the published literature! Unfortunately it sounds like the mistake is too small to matter much... but still, it's amusing. What next? Is there anything I can do to help? I could help polish the page [[Another look at climate sensitivity]], but I don't understand the graphs at the bottom: which graph came from which procedure. If someone writes a quick note to make that clear, I can add more detail.`

I'm not bummed about finding a mistake. It happens. I'm just a bit bummed that the mistake is significant enough to void much of the discussion in the paper.

I admit that I have some preconceived notions that climate science was kind of soft to begin with and easy-to-find-mistakes like that enforce the notion and imply the publishing standards are pretty poor.

`I'm not bummed about finding a mistake. It happens. I'm just a bit bummed that the mistake is significant enough to void much of the discussion in the paper. I admit that I have some preconceived notions that climate science was kind of soft to begin with and easy-to-find-mistakes like that enforce the notion and imply the publishing standards are pretty poor.`

Eric said:

I wouldn't be surprized if this kind of mistake could get passed the referees in other fields too.

`Eric said: > imply the publishing standards are pretty poor I wouldn't be surprized if this kind of mistake could get passed the referees in other fields too.`

This is an old issue about what referees are supposed to do (at least in experimental sciences, in "pure logic" sciences it might be different). My view, and the instructions I received back when I did refereeing, was that refereeing is supposed to determine if the experimental procedure seems sound, if the investigation is novel and "significant" (which is the most moveable feast criterion), and if the results that were claimed to be obtained have been interpreted appropriately. With the exception of papers that were clearly so poorly put together that the likelihood of results that do not accurately reflect what would happen if you followed the experimental description (whether from typos/mistakes or something else) looked appreciable, it's not expect that (under the current system) the referee will confirm that the results are indeed obtained from the experiment. That is supposed to come from independent replication of other groups on their way to trying to build on these results.

Of course, computer simulation papers are on the borderline where it's just about possible to repeat the experiments in some cases, such as this one, and maybe different standards should be settled upon for those cases.

`This is an old issue about what referees are supposed to do (at least in experimental sciences, in "pure logic" sciences it might be different). My view, and the instructions I received back when I did refereeing, was that refereeing is supposed to determine if the experimental procedure seems sound, if the investigation is novel and "significant" (which is the most moveable feast criterion), and if the results that were claimed to be obtained have been interpreted appropriately. With the exception of papers that were clearly so poorly put together that the likelihood of results that do not accurately reflect what would happen if you followed the experimental description (whether from typos/mistakes or something else) looked appreciable, it's not expect that (under the current system) the referee will confirm that the results are indeed obtained from the experiment. That is supposed to come from independent replication of other groups on their way to trying to build on these results. Of course, computer simulation papers are on the borderline where it's just about possible to repeat the experiments in some cases, such as this one, and maybe different standards should be settled upon for those cases.`

Eric, why do you say the mistake voids much of their discussion? The qualitative structure of the bifurcation doesn't change. The value $\mu = 1$, which they want to correspond to our current climate, leads in the corrected plot only to a cold-Earth scenario, so it's clear that these parameter choices don't fit the real world very well. But I'm not sure how much this should bother me.

I think understanding the paper really requires putting it in a broader context, which I haven't had the time to do the reading for yet. Being new here, I don't want to derail the ongoing discussion, which I have only a fuzzy sense of the goal of. But the questions that come to mind when I read the paper are:

How reliable are the toy models of ice-albedo feedback and the greenhouse effect? Should I believe they accurately capture at least semi-quantitatively the dynamics of these mechanisms, so that plugging them into a simple energy balance model is not a crazy model for the state of the Earth?

If these are good toy models, to what extent are their parameters known? Can we pick different, reasonable values that produce a new bifurcation diagram that passes near a real-world temperature value on the "warm" branch, and exhibits the same qualitative structure?

To what extent should I believe that the existence of bifurcations in a one-dimensional model (i.e. an energy balance model that averages over the whole Earth) is a good guide to the dynamics of the whole system? It suggests pretty dramatic effects, like hysteresis, that I'm not sure I believe would apply to a more complete climate model. It might be that a more complicated model, for instance, would lead to dynamics where the globe separates into regions of different temperatures in a way that makes the averaged curve monotonic. (I'm thinking by loose analogy to how effective potentials in physics are always convex, even in the presence of phase transitions, because minimum-energy states don't have to be homogeneous but can have regions of different phases. If this analogy doesn't help anyone else, please ignore it.)

The part of the paper that I was bummed about was the beginning, where they discuss Roe & Baker's paper on feedbacks and climate sensitivity. The Roe & Baker paper seems to have gotten a lot of attention; I was aware of it because RealClimate blogged about it when it came out. If Zaliapin and Ghil's discussion of it represents it accurately, this is pretty depressing, since it sounds like Roe and Baker made trivial mistakes in thinking about what Taylor series mean! Of course, Roe and Baker have written a Comment in response, and Zaliapin and Ghil have written a Reply to the comment. Sorting out who is correct on this basic point about series expansions and their implication for uncertainties in outputs with feedback versus uncertainties in inputs should be easy to do (though it might turn out that their real disagreement isn't about trivial mathematics at all, but rather about reasonable ranges for the values of various parameters).

`Eric, why do you say the mistake voids much of their discussion? The qualitative structure of the bifurcation doesn't change. The value $\mu = 1$, which they want to correspond to our current climate, leads in the corrected plot only to a cold-Earth scenario, so it's clear that these parameter choices don't fit the real world very well. But I'm not sure how much this should bother me. I think understanding the paper really requires putting it in a broader context, which I haven't had the time to do the reading for yet. Being new here, I don't want to derail the ongoing discussion, which I have only a fuzzy sense of the goal of. But the questions that come to mind when I read the paper are: 1. How reliable are the toy models of ice-albedo feedback and the greenhouse effect? Should I believe they accurately capture at least semi-quantitatively the dynamics of these mechanisms, so that plugging them into a simple energy balance model is not a crazy model for the state of the Earth? 2. If these are good toy models, to what extent are their parameters known? Can we pick different, reasonable values that produce a new bifurcation diagram that passes near a real-world temperature value on the "warm" branch, and exhibits the same qualitative structure? 3. To what extent should I believe that the existence of bifurcations in a one-dimensional model (i.e. an energy balance model that averages over the whole Earth) is a good guide to the dynamics of the whole system? It suggests pretty dramatic effects, like hysteresis, that I'm not sure I believe would apply to a more complete climate model. It might be that a more complicated model, for instance, would lead to dynamics where the globe separates into regions of different temperatures in a way that makes the averaged curve monotonic. (I'm thinking by loose analogy to how effective potentials in physics are always convex, even in the presence of phase transitions, because minimum-energy states don't have to be homogeneous but can have regions of different phases. If this analogy doesn't help anyone else, please ignore it.) The part of the paper that I was bummed about was the beginning, where they discuss Roe & Baker's paper on feedbacks and climate sensitivity. The Roe & Baker paper seems to have gotten a lot of attention; I was aware of it because RealClimate blogged about it when it came out. If Zaliapin and Ghil's discussion of it represents it accurately, this is pretty depressing, since it sounds like Roe and Baker made trivial mistakes in thinking about what Taylor series mean! Of course, Roe and Baker have written a [Comment in response](http://www.nonlin-processes-geophys.net/18/125/2011/npg-18-125-2011.pdf), and Zaliapin and Ghil have written a [Reply to the comment](http://www.nonlin-processes-geophys.net/18/129/2011/npg-18-129-2011.pdf). Sorting out who is correct on this basic point about series expansions and their implication for uncertainties in outputs with feedback versus uncertainties in inputs should be easy to do (though it might turn out that their real disagreement isn't about trivial mathematics at all, but rather about reasonable ranges for the values of various parameters).`

Sorry for the curt responses. I'm not very emotional about this paper, just very short of time. It's no big deal. Basically, I wasn't very enthusiastic about the paper when I first saw it, but decided to work through it anyway hoping to learn something. I don't think I learned much and the paper contains an irritating error that demonstrates a frustrating carelessness on the part of the author.

Anyway, I think there are cool things to be done and I still hope to be productive. Best wishes.

`Sorry for the curt responses. I'm not very emotional about this paper, just very short of time. It's no big deal. Basically, I wasn't very enthusiastic about the paper when I first saw it, but decided to work through it anyway hoping to learn something. I don't think I learned much and the paper contains an irritating error that demonstrates a frustrating carelessness on the part of the author. Anyway, I think there are cool things to be done and I still hope to be productive. Best wishes.`

David said:

I think it would be significant progress if the authors published and versioned their code along with the paper, so that at least they themselves can reproduce what they have done, and others are able to confirm their suspicians by looking at the published code.

Matt said:

The overall goal is to gain a better understanding of the most simple models used in climate science, by reproducing the results of published papers - this is a first step into more serious investigations. Part of this is reproducing computer generated graphics.

I have a hidden agenda here, which is:

a) demonstrate how the software part is done with modern tools,

b) demonstrate how the code should be structured and published along with results.

But that's just me, the overall goal is pretty much thinking and learning about the questions you asked.

Climate science is very demanding, because it is highly inter-disciplinary. That means that a lot of people will make a lot of mistakes that look pretty darn stupid to the experts, when they write about climate science. This does not depress me, I see it as a call to duty.

Physicists and mathematicians are no exception, by the way, I've seen peer-reviewed papers published in physics journals by tenured physicists from major German universisties that are - well - embarrasing. The notorious Gerlich-Tscheuschner paper comes to my mind.

`David said: <blockquote> <p> Of course, computer simulation papers are on the borderline where it's just about possible to repeat the experiments in some cases, such as this one, and maybe different standards should be settled upon for those cases. </p> </blockquote> I think it would be significant progress if the authors published and versioned their code along with the paper, so that at least they themselves can reproduce what they have done, and others are able to confirm their suspicians by looking at the published code. Matt said: <blockquote> <p> Being new here, I don't want to derail the ongoing discussion, which I have only a fuzzy sense of the goal of. </p> </blockquote> The overall goal is to gain a better understanding of the most simple models used in climate science, by reproducing the results of published papers - this is a first step into more serious investigations. Part of this is reproducing computer generated graphics. I have a hidden agenda here, which is: a) demonstrate how the software part is done with modern tools, b) demonstrate how the code should be structured and published along with results. But that's just me, the overall goal is pretty much thinking and learning about the questions you asked. <blockquote> <p> ...this is pretty depressing, since it sounds like Roe and Baker made trivial mistakes in thinking about what Taylor series mean! </p> </blockquote> Climate science is very demanding, because it is highly inter-disciplinary. That means that a lot of people will make a lot of mistakes that look pretty darn stupid to the experts, when they write about climate science. This does not depress me, I see it as a call to duty. Physicists and mathematicians are no exception, by the way, I've seen peer-reviewed papers published in physics journals by tenured physicists from major German universisties that are - well - embarrasing. The notorious <a href="http://arxiv.org/abs/0707.1161">Gerlich-Tscheuschner</a> paper comes to my mind.`

I just added a summary table of tipping points on Tipping point from Thompson and Sieber excellent survey. It needed to get some perspective on which tipping points /comb of bifurcations.

Yes I was a bit bummed but i implemented most of the things we talked about , and also thought - like Matt - that compared to what ZG10 where bashing it was a minor thing. I think that ZG is correct about the Taylor expansion, but I leave that you mathematicians to resolve

So my "bumminess" had more todo that I needed to get some fresh input and get the discussion going here for details on how to proceed. As the last part i the Sage worksheet is simply an interactive feature - you cannot see that but I will take some static snapshots so you can see.

`I just added a summary table of tipping points on [[Tipping point]] from Thompson and Sieber excellent survey. It needed to get some perspective on which tipping points /comb of bifurcations. Yes I was a bit bummed but i [implemented most of the things we talked about](http://demo.sagenb.org/home/pub/92/) , and also thought - like Matt - that compared to what ZG10 where bashing it was a minor thing. I think that ZG is correct about the Taylor expansion, but I leave that you mathematicians to resolve So my "bumminess" had more todo that I needed to get some fresh input and get the discussion going here for details on how to proceed. As the last part i the Sage worksheet is simply an interactive feature - you cannot see that but I will take some static snapshots so you can see.`

John Baez said

It is that the scientific spirit? Or an academic careerist spirit? I know you are not an academic careerist. However, there is a fine line between being pleased that you found a mistake, and being pleased that someone else made a mistake that you could find. In my opinion you should

notbe happy that they are wrong, though it of course nice for Azimuth that people here were the ones to find it.`John Baez said > All this sounds great except 4). You should be _happy_ because they're wrong - and you caught it! Where's your scientific spirit? There's nothing more fun than finding a mistake in someone else's work! I _love_ it when that happens! It is that the scientific spirit? Or an academic careerist spirit? I know you are not an academic careerist. However, there is a fine line between being pleased that you found a mistake, and being pleased that someone else made a mistake that you could find. In my opinion you should _not_ be happy that they are wrong, though it of course nice for Azimuth that people here were the ones to find it.`

David said:

and Tim replied:

I agree. Some journals are beginning to insist on this. One thing I like about R is its support for 'reproducible research'. From CRAN = The Comprehensive R Archive Network:

Another link: Reproducible Research with R, LATEX, & Sweave

`David said: > Of course, computer simulation papers are on the borderline where it's just about possible to repeat the experiments in some cases, such as this one, and maybe different standards should be settled upon for those cases. and Tim replied: > I think it would be significant progress if the authors published and versioned their code along with the paper, so that at least they themselves can reproduce what they have done, and others are able to confirm their suspicians by looking at the published code. I agree. Some journals are beginning to insist on this. One thing I like about R is its support for 'reproducible research'. From CRAN = The Comprehensive R Archive Network: > The goal of [reproducible research](http://cran.r-project.org/web/views/ReproducibleResearch.html) is to tie specific instructions to data analysis and experimental data so that scholarship can be recreated, better understood and verified. Another link: [Reproducible Research with R, LATEX, & Sweave](http://biostat.mc.vanderbilt.edu/wiki/pub/Main/TheresaScott/ReproducibleResearch.TAScott.handout.pdf)`

These are certainly advances, but it's worth bearing in mind that actively debugging (going looking for non-obvious problems) a large codebase is an intensive job. (Suppose I'm using a version of BLAS that generates incorrect results, or the compiler is miscompiling something because it's doing a wrong analysis, etc.)

`These are certainly advances, but it's worth bearing in mind that actively debugging (going looking for non-obvious problems) a large codebase is an intensive job. (Suppose I'm using a version of BLAS that generates incorrect results, or the compiler is miscompiling something because it's doing a wrong analysis, etc.)`

It certainly is, but more often than compilers or BLAS going wrong, somebody just makes a sign error.

`> actively debugging (going looking for non-obvious problems) [in] a large codebase is an intensive job. It certainly is, but more often than compilers or BLAS going wrong, somebody just makes a sign error.`

Sorry again for falling into negative territory. Getting back to the positives:

John wrote:

Although I presume you are not much into writing code, I think it would be nice if you could at least run any code we produce here. Seeing the code in action might help generate ideas.

Since I already made all the obvious mistakes, I can help you avoid them. I suggest you install the following bare-bones stuff:

With the direction I see you going in your current research, at some point I think you will reach the limits of what you can do without computation anyway so it is probably a good idea you start getting your hands dirty now.

`Sorry again for falling into negative territory. Getting back to the positives: John wrote: >What next? Is there anything I can do to help? Although I presume you are not much into writing code, I think it would be nice if you could at least run any code we produce here. Seeing the code in action might help generate ideas. Since I already made all the obvious mistakes, I can help you avoid them. I suggest you install the following bare-bones stuff: * [Python](http://www.python.org/getit/) - I'd go for 2.7 for now. Upgrading to 3.x later should be simple, but some add-ons are not compatible with 3.x yet. * [Numpy and Scipy](http://www.scipy.org/Download) - Make sure you download the version corresponding to the Python version you downloaded. * [matplotlib](http://sourceforge.net/projects/matplotlib/files/matplotlib/) - Again, make sure the version you download corresponds to the Python version. With the direction I see you going in your current research, at some point I think you will reach the limits of what you can do without computation anyway so it is probably a good idea you start getting your hands dirty now.`

Instead of trying to worry about which is the right version of all of those, you can take the lazy route and get the Enthought python distribution, which puts it all in one convenient package (free for individual academic use). It's worked well for me, at least.

`Instead of trying to worry about which is the right version of all of those, you can take the lazy route and get the [Enthought python distribution](http://www.enthought.com/products/epd.php), which puts it all in one convenient package (free for individual academic use). It's worked well for me, at least.`

True. Enthought looks like a good option for academic purposes. Unfortunately, I'm not associated with a university. I almost wish we could keep everything purely open. I think others here are also not affiliated with a university.

`True. Enthought looks like a good option for academic purposes. Unfortunately, I'm not associated with a university. I almost wish we could keep everything purely open. I think others here are also not affiliated with a university.`

I'm not affiliated with a university, but besides that, I believe that implementations of algorithms commonly used in statistical analysis etc. should be free and open source, and there are enough people out there who think likewise and already have put a lot of work into open source, free implementations. Our best example at this point is the Sage project.

Of course I've used Maple for my diploma thesis, to solve a large system of linear equations in closed form: If you have limited time to produce a certain result and some commercial software provided by your university is the best choice, I won't blame you if you use it :-)

`I'm not affiliated with a university, but besides that, I believe that implementations of algorithms commonly used in statistical analysis etc. should be free and open source, and there are enough people out there who think likewise and already have put a lot of work into open source, free implementations. Our best example at this point is the Sage project. Of course I've used Maple for my diploma thesis, to solve a large system of linear equations in closed form: If you have limited time to produce a certain result and some commercial software provided by your university is the best choice, I won't blame you if you use it :-)`

Eric wrote:

That would be

— it would mean we could easily publish a paper correcting theirs! But alas, I don't see that this mistake voidsgreatanythingin their paper. It changes the shape of one graph slightly, but that's all. And all that matters about this graph is its general shape: catastrophe theory is robust against small perturbations. So nothing significant changes.Am I confused? Maybe I missed some other consequence of their mistake?

`Eric wrote: > I'm just a bit bummed that the mistake is significant enough to void much of the discussion in the paper. That would be _**great**_ — it would mean we could easily publish a paper correcting theirs! But alas, I don't see that this mistake voids _anything_ in their paper. It changes the shape of one graph slightly, but that's all. And all that matters about this graph is its general shape: catastrophe theory is robust against small perturbations. So nothing significant changes. Am I confused? Maybe I missed some other consequence of their mistake?`

John Baez said:

Graham Jones wrote:

I was being ironic, of course. It's just very amusing for me to hear people being "bummed" by finding a mistake in someone else's paper, since most academics I know are quite pleased to catch errors in other people's work.

It's true that there's a fine line here. But there actually is something beneficial about having lots of ambitious careerists eager to find and catch errors: it means errors tend to get caught!

`John Baez said: > All this sounds great except 4). You should be happy because they're wrong - and you caught it! Where's your scientific spirit? There's nothing more fun than finding a mistake in someone else's work! I love it when that happens! Graham Jones wrote: > It is that the scientific spirit? Or an academic careerist spirit? I know you are not an academic careerist. However, there is a fine line between being pleased that you found a mistake, and being pleased that someone else made a mistake that you could find. In my opinion you should not be happy that they are wrong, though it of course nice for Azimuth that people here were the ones to find it. I was being ironic, of course. It's just very amusing for me to hear people being "bummed" by finding a mistake in someone else's paper, since most academics I know are quite pleased to catch errors in other people's work. It's true that there's a fine line here. But there actually is something beneficial about having lots of ambitious careerists eager to find and catch errors: it means errors tend to get caught!`

Ok. Back in this comment, I wrote:

Here are the results...

I started with a tiny amplitude $a = .01$:

The little red blurb in the third subplot denotes the temperature variation around equilibrium.

Then I bumped up the amplitude a bit to $a = .05$:

Lo and behold, the temperature suddenly jumps down 120 degrees and remains stable there.

The picture looks similar with $a = .1$:

Although a bit extreme, here is $a = .5$:

It looks like it doesn't take much to launch us into an ice age and it seem next to impossible to warm up once we're in an ice age. Or so it seems to this layman...

Here is the code:

I'm new to Python and I'm sure this is a total hack job, so if anyone has any tips for writing this more elegantly, please let me know.

There could be bugs as well, so comments welcome.

`Ok. Back in [this comment](http://www.math.ntnu.no/~stacey/Mathforge/Azimuth/comments.php?DiscussionID=569&Focus=3621#Comment_3621), I wrote: >An idea I have to take this in the desired direction of a dynamic insolation driver is to >1. Get the steady state correct (which means getting this chart right). >1. Begin at some point on the chart that is in steady state. This will be our initial conditions. >1. Turn on a driver $\mu(t) = \mu_0 + u(t) a sin(\omega t)$, where $u(t)$ is the unit step function. >1. Solve the time evolution of $T$ numerically. Here are the results... I started with a tiny amplitude $a = .01$: <img src="http://www.azimuthproject.org/azimuth/files/ZG_FD_01.png" width = "500" alt = ""/> The little red blurb in the third subplot denotes the temperature variation around equilibrium. Then I bumped up the amplitude a bit to $a = .05$: <img src="http://www.azimuthproject.org/azimuth/files/ZG_FD_05.png" width = "500" alt = ""/> Lo and behold, the temperature suddenly jumps down 120 degrees and remains stable there. The picture looks similar with $a = .1$: <img src="http://www.azimuthproject.org/azimuth/files/ZG_FD_10.png" width = "500" alt = ""/> Although a bit extreme, here is $a = .5$: <img src="http://www.azimuthproject.org/azimuth/files/ZG_FD_50.png" width = "500" alt = ""/> It looks like it doesn't take much to launch us into an ice age and it seem next to impossible to warm up once we're in an ice age. Or so it seems to this layman... Here is the code: import numpy import scipy import matplotlib from pylab import * sigma = 5.6697e-8; Q0 = 342.5; kappa= 0.05; def alpha(kappa,T): c1 = 0.15; c2 = 0.7; Tc = 273.0; return c1 + .5*c2*(1 - tanh(kappa*(T - Tc))) def g(T): m = 0.4; T0 = (1.9e-15)**(-1/6.); return 1 - m*tanh((T/T0)**6) def Ro(sigma,T): return sigma*g(T)*T**4 def Ri(mu,Q0,kappa,T): return mu*Q0*(1 - alpha(kappa,T)) def mu_eq(sigma,Q0,kappa,T): return Ro(sigma,T)/Ri(1.,Q0,kappa,T) def get_mu(mu0,a,omega,t): mu = 1.*t; mu[t < 0] = mu0; mu[t >= 0] = mu0+a*sin(omega*t[t >=0]) return mu T_init = 300.; mu_init = mu_eq(sigma,Q0,kappa,T_init); a = .5; omega = .1; c = 1.; nt = 100.; dt = 2*pi/(nt*omega); nmin = -nt; nmax = 9*nt; time = arange(nmin, nmax, dt); mu = get_mu(mu_init,a,omega,time); Rnet = 1.*time; T = 1.*time; i = 0; T_prev = 1.*T_init; for t in time: Rnet[i] = Ri(mu[i],Q0,kappa,T_prev)-Ro(sigma,T_prev) T[i] = T_prev+(dt/c)*Rnet[i] T_prev = T[i] i = i+1; subplot(311) plot(time, mu,'b-', lw=2) grid(True) xlabel('Time') ylabel('mu') subplot(312) plot(time, T,'b-', lw=2) grid(True) xlabel('Time') ylabel('Temperature') subplot(313) Tp = arange(100.0, 400.0+1.0, 1.0) plot(Tp, mu_eq(sigma,Q0,kappa,Tp), 'b-', T, mu, 'r-', lw=2) grid(True) xlabel('Temperature') ylabel('mu') show() I'm new to Python and I'm sure this is a total hack job, so if anyone has any tips for writing this more elegantly, please let me know. There could be bugs as well, so comments welcome.`

Hi Eric, it might be an idea to do the development on a wiki-page, with some pointers comments here, so that the work is more permanently recorded.

`Hi Eric, it might be an idea to do the development on a wiki-page, with some pointers comments here, so that the work is more permanently recorded.`