I mentioned that I would place a canonical representation of the SOI model on the forum. This is a typical combination of the factors representing the QBO, TSI, Chandler Wobble (CW), plus a slowly varying multi-decadal term (MW). Each of these is represented as a RHS forcing and as a Mathieu/Hill modulation term. The characteristic frequency of the wave equation is embodied in the CF term.
(* The SOIM differential equation solved with the Mathematica NDSolve function *)
RHS[x_]:= .9*(.063*qbo[x+.51]*(1+.108*tsi[x-1.1])-.0192*cw[x+1.7]-.012*tsi[x-4.2])+.018*mw[x];
Hill[x_] := 1.32*(.37*tsi[x-1.6]-.52*cw[x+2.94]-.24*qbo[x+.6]*(1-.54*tsi[x-3.1])) ;
SOLN = NDSolve[{y''[x]+(CF+Hill[x])*y[x] == RHS[x], y'[52]==.0, y[52]==.0015}, y, {x, -20, 160}];
That's it, with initial conditions specified and x in years.
These are the fundamental frequencies for each of the factors. The characteristic frequency was determined by Clarke at FSU.
(* constants for QBO=2.33yr, TS=10.6yr, CW=6.46yr(432day beat with 365day) and the Clarke CharFreq=4.24yr *)
QB=2.696; TS:=0.5909; CF = 2.197; CW = 0.9726; MW := 0.116;
For the terms representing the four factors, I did my best to emulate the real data as a combination of sinusoidal signals leading to the following closed-form expressions. This makes it easier to solve the differential equation -- as the Mathematica solver can interpolate as necessary to converge to a solution:
(* closed form representations of the QBO, TSI, CW, MW data *)
qbo[x_] := 1.1*Cos[QB*x+1.5+.35*Sin[.515*x-2.] ]; (* QBO time series *)
tsi[x_] := (Cos[TS*x+.357]*(1+3.5*Sin[(x-35.4)/217]) + (* TSI time series low-pass filtered *)
.6*sig2[x,5.7,25.1,28.]+ 2.5*sig2[x,1.6,87.9,89.6])*((1+.3*sigmoid[x,106.4])/1.3);
cw[x_] := 1.11*Sin[CW*(x+0.88)+0.2*Sin[0.361*(x+5.3)]] ; (* Chandler wobble time series *)
mw[x_] := Cos[MW*x-0.7]; (* inferred multi-decadal forcing, LOD ? Markowitz wobble ? *)
The following figure shows how well the waveforms match the actual data for QBO, CW, and TSI. The correlation coefficient multiplied by 100 is shown above each time-series. The closer one can correlate to these waveforms , the better the fit to the SOI data.
And here is the fit of the canonical model to the SOI data for the years ranging from 1880 to 2013:
The bottom two graphs are the wavelet scalograms of the data and of the model. I am still looking for a goodness of fit measure for these scalograms.
The bar chart is the classifier result. Calculated for every month of data, it shows the probability that the model and SOI data have the same sign, which is above 70% of the time. Although this may not look so good on first glance, it is excellent considering that the two dipole values of Tahiti and Darwin only have the opposite signed excursion only about 2/3 of the time! That is what we are dealing with in terms of measuring the actual climate dipole in action.
---
*EDIT*
Jim reminded me that the sigmoid utility functions were inadvertently omitted
(* shaping functions for capturing model of TSI *)
sig2[x_, W_, T1_, T2_] := -1/(1+Exp[W*(x-T1)])+1/(1+Exp[W*(x-T2)]);
sigmoid[x_, T_] := 1/(1+Exp[2.7*(x-T)]);
(* The SOIM differential equation solved with the Mathematica NDSolve function *)
RHS[x_]:= .9*(.063*qbo[x+.51]*(1+.108*tsi[x-1.1])-.0192*cw[x+1.7]-.012*tsi[x-4.2])+.018*mw[x];
Hill[x_] := 1.32*(.37*tsi[x-1.6]-.52*cw[x+2.94]-.24*qbo[x+.6]*(1-.54*tsi[x-3.1])) ;
SOLN = NDSolve[{y''[x]+(CF+Hill[x])*y[x] == RHS[x], y'[52]==.0, y[52]==.0015}, y, {x, -20, 160}];
That's it, with initial conditions specified and x in years.
These are the fundamental frequencies for each of the factors. The characteristic frequency was determined by Clarke at FSU.
(* constants for QBO=2.33yr, TS=10.6yr, CW=6.46yr(432day beat with 365day) and the Clarke CharFreq=4.24yr *)
QB=2.696; TS:=0.5909; CF = 2.197; CW = 0.9726; MW := 0.116;
For the terms representing the four factors, I did my best to emulate the real data as a combination of sinusoidal signals leading to the following closed-form expressions. This makes it easier to solve the differential equation -- as the Mathematica solver can interpolate as necessary to converge to a solution:
(* closed form representations of the QBO, TSI, CW, MW data *)
qbo[x_] := 1.1*Cos[QB*x+1.5+.35*Sin[.515*x-2.] ]; (* QBO time series *)
tsi[x_] := (Cos[TS*x+.357]*(1+3.5*Sin[(x-35.4)/217]) + (* TSI time series low-pass filtered *)
.6*sig2[x,5.7,25.1,28.]+ 2.5*sig2[x,1.6,87.9,89.6])*((1+.3*sigmoid[x,106.4])/1.3);
cw[x_] := 1.11*Sin[CW*(x+0.88)+0.2*Sin[0.361*(x+5.3)]] ; (* Chandler wobble time series *)
mw[x_] := Cos[MW*x-0.7]; (* inferred multi-decadal forcing, LOD ? Markowitz wobble ? *)
The following figure shows how well the waveforms match the actual data for QBO, CW, and TSI. The correlation coefficient multiplied by 100 is shown above each time-series. The closer one can correlate to these waveforms , the better the fit to the SOI data.
And here is the fit of the canonical model to the SOI data for the years ranging from 1880 to 2013:
The bottom two graphs are the wavelet scalograms of the data and of the model. I am still looking for a goodness of fit measure for these scalograms.
The bar chart is the classifier result. Calculated for every month of data, it shows the probability that the model and SOI data have the same sign, which is above 70% of the time. Although this may not look so good on first glance, it is excellent considering that the two dipole values of Tahiti and Darwin only have the opposite signed excursion only about 2/3 of the time! That is what we are dealing with in terms of measuring the actual climate dipole in action.
---
*EDIT*
Jim reminded me that the sigmoid utility functions were inadvertently omitted
(* shaping functions for capturing model of TSI *)
sig2[x_, W_, T1_, T2_] := -1/(1+Exp[W*(x-T1)])+1/(1+Exp[W*(x-T2)]);
sigmoid[x_, T_] := 1/(1+Exp[2.7*(x-T)]);
Version 2.1.8p2
