Application Center - Maplesoft
  Submit your work
Application Center Applications Global Temperature Anomaly Preview

App Preview:

Global Temperature Anomaly

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



> 

 

 

January  2016

 

                                           Global Temperature Anomaly

 

                                         Univ.-Prof. Dr.-Ing. habil. Josef  BETTEN

                                                   RWTH Aachen University

                        Mathematical Models in Materials Science and Continuum Mechanics

                                                           Augustinerbach 4-20

                                                      D-52056  A a c h e n ,  Germany

 

                                                      <betten@mmw.rwth-aachen.de>

 

Abstract

 

The temperature anomaly or temperature index is defined as the change from a reference temperature or a long-term mean value. A positive ( negative) anomaly indicates that a measured temperuture is warmer (cooler) than the reference value. In this worksheet anomalies have been based upon the period between 1951 to 1980. We consider especially  temperature change from 1980 to 2015. Since the global temperatur is increasing in the last years, i.e., the temperature anomaly is positive in the following. One can differentiate between global land-ocean or surface temperature index ( °C ), for instance.

 

Global Land-Ocean Temperature Index  ( °C )

 

The following DATA are based on experiments do to Nasa   (Anomaly with base: 1951 - 1980). These experimental values have been interpolated  by a  cubic spline, which is approximated by  FOURIER-series.  

 

> 

restart:
> 

with(stats):
> 

with(CurveFitting):
> 

DATA:=               [1980,0.28],[1981,0.33],[1982,0.13],[1983,0.31],[1984,0.16], [1985,0.12],[1986,0.19],[1987,0.34],[1988,0.40],[1989,0.29], [1990,0.44],[1991,0.43],[1992,0.23],[1993,0.24],[1994,0.32], [1995,0.46],[1996,0.35],[1997,0.48],[1998,0.64],[1999,0.42], [2000,0.42],[2001,0.55],[2002,0.63],[2003,0.62],[2004,0.55], [2005,0.69],[2006,0.63],[2007,0.66],[2008,0.54],[2009,0.65], [2010,0.72],[2011,0.61],[2012,0.63],[2013,0.66],[2014,0.75],          [2015,0.85];

DATA := [1980, .28], [1981, .33], [1982, .13], [1983, .31], [1984, .16], [1985, .12], [1986, .19], [1987, .34], [1988, .40], [1989, .29], [1990, .44], [1991, .43], [1992, .23], [1993, .24], [1994, .32], [1995, .46], [1996, .35], [1997, .48], [1998, .64], [1999, .42], [2000, .42], [2001, .55], [2002, .63], [2003, .62], [2004, .55], [2005, .69], [2006, .63], [2007, .66], [2008, .54], [2009, .65], [2010, .72], [2011, .61], [2012, .63], [2013, .66], [2014, .75], [2015, .85]

 (1)
> 

# Cubic Spline Function:
> 

Sp(x):=Spline([DATA],x,degree=3):
> 

# This function is not printed because of its length. One can print it, if we insert a semicolon (;) instead of the doppel point (:) at the end of the spline  comand. The following Figure shows the result of the Cubic Spline Interpolation.

> 

alias(th=thickness,co=color):
> 

p[1]:=plot([DATA],x=1980..2015,th=3,co=black, axes=boxed,style=point,symbol=cross,symbolsize=40):
> 

p[2]:=plot(Sp(x),x=1980..2015,th=3,co=black,                        title="Global Temperature Index (°C)"):
> 

p[3]:=plots[textplot]({[1990,0.75,`DATA by NASA`],              [2007,0.25,`Cubic Spline Interpolation`]}):
> 

plots[display]({seq(p[k],k=1..3)});

 
> 

 

 

Approximation of the Splinefunction by FOURIER-Series

 

First let us map linear the annual abscissa  X = [1980...2015]  to the period  x = [0...2*Pi] according to:

> 

restart:

> 

x:=2*Pi*(X-1980)/35;

x := 2*Pi*(X-1980)/35

 (2)
> 

# The experimental data are  expressed as:
> 

for i from 1980 to 2015 do x[i]:=[data]=2*Pi*(i-1980)/35 od:
> 

datax:=array([seq(2*Pi*(i-1980)/35,i=1980..2015)]);

datax := Vector[row](36, {(1) = 0, (2) = (2/35)*Pi, (3) = (4/35)*Pi, (4) = (6/35)*Pi, (5) = (8/35)*Pi, (6) = (2/7)*Pi, (7) = (12/35)*Pi, (8) = (2/5)*Pi, (9) = (16/35)*Pi, (10) = (18/35)*Pi, (11) = (4/7)*Pi, (12) = (22/35)*Pi, (13) = (24/35)*Pi, (14) = (26/35)*Pi, (15) = (4/5)*Pi, (16) = (6/7)*Pi, (17) = (32/35)*Pi, (18) = (34/35)*Pi, (19) = (36/35)*Pi, (20) = (38/35)*Pi, (21) = (8/7)*Pi, (22) = (6/5)*Pi, (23) = (44/35)*Pi, (24) = (46/35)*Pi, (25) = (48/35)*Pi, (26) = (10/7)*Pi, (27) = (52/35)*Pi, (28) = (54/35)*Pi, (29) = (8/5)*Pi, (30) = (58/35)*Pi, (31) = (12/7)*Pi, (32) = (62/35)*Pi, (33) = (64/35)*Pi, (34) = (66/35)*Pi, (35) = (68/35)*Pi, (36) = 2*Pi})

 (3)
> 

datax:=evalf(array([seq(2*Pi*(i-1980)/35,i=1980..2015)]),3);

datax := Vector[row](36, {(1) = 0., (2) = .179, (3) = .358, (4) = .537, (5) = .719, (6) = .898, (7) = 1.08, (8) = 1.26, (9) = 1.43, (10) = 1.61, (11) = 1.79, (12) = 1.98, (13) = 2.15, (14) = 2.33, (15) = 2.51, (16) = 2.69, (17) = 2.87, (18) = 3.05, (19) = 3.23, (20) = 3.42, (21) = 3.58, (22) = 3.77, (23) = 3.96, (24) = 4.11, (25) = 4.30, (26) = 4.49, (27) = 4.68, (28) = 4.84, (29) = 5.02, (30) = 5.21, (31) = 5.37, (32) = 5.56, (33) = 5.75, (34) = 5.93, (35) = 6.09, (36) = 6.28})

 (4)
> 

whattype(datax);

symbol

 (5)
> 

datay:=array([0.28,0.33,0.13,0.31,0.16,0.12,0.19,0.34,                 0.40,0.29, 0.44,0.43,0.23,0.24,0.32,0.46,0.35,0.48,                0.64,0.42,0.42,0.55, 0.63,0.62,0.55,0.69,0.63,0.66,     0.54,0.65,0.72,0.61,0.63,0.66,0.75,0.85]);

datay := Vector[row](36, {(1) = .28, (2) = .33, (3) = .13, (4) = .31, (5) = .16, (6) = .12, (7) = .19, (8) = .34, (9) = .40, (10) = .29, (11) = .44, (12) = .43, (13) = .23, (14) = .24, (15) = .32, (16) = .46, (17) = .35, (18) = .48, (19) = .64, (20) = .42, (21) = .42, (22) = .55, (23) = .63, (24) = .62, (25) = .55, (26) = .69, (27) = .63, (28) = .66, (29) = .54, (30) = .65, (31) = .72, (32) = .61, (33) = .63, (34) = .66, (35) = .75, (36) = .85})

 (6)
> 

whattype(datay);

symbol

 (7)
> 

DATA:=seq([datax[i],datay[i]],i=1..36);

DATA := [0., .28], [.179, .33], [.358, .13], [.537, .31], [.719, .16], [.898, .12], [1.08, .19], [1.26, .34], [1.43, .40], [1.61, .29], [1.79, .44], [1.98, .43], [2.15, .23], [2.33, .24], [2.51, .32], [2.69, .46], [2.87, .35], [3.05, .48], [3.23, .64], [3.42, .42], [3.58, .42], [3.77, .55], [3.96, .63], [4.11, .62], [4.30, .55], [4.49, .69], [4.68, .63], [4.84, .66], [5.02, .54], [5.21, .65], [5.37, .72], [5.56, .61], [5.75, .63], [5.93, .66], [6.09, .75], [6.28, .85]

 (8)
> 

with(stats):   with(CurveFitting):
> 

# The cubic splinefunction is given by:
> 

Sp(x):=Spline([DATA],x,degree=3):
> 

# The results are represented in the following Figure:
> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=40,                             title="Global Temperature Index (°C)"):
> 

p[3]:=plots[textplot]({[2,0.75,`DATA by Nasa`],                      [4.5,0.25,`Cubic Spline Interpolation`]}):
> 

plots[display](seq(p[k],k=1..3));

 
> 

 

> 

FOURIER_series(x):=          a[0]/2+sum(a[k]*cos(k*x)+b[k]*sin(k*x),k=1..infinity);

FOURIER_series(x) := (1/2)*a[0]+sum(a[k]*cos(k*x)+b[k]*sin(k*x), k = 1 .. infinity)

 (9)
> 

a[k]:=(1/Pi)*Int(f(x)*cos(k*x),x=0..2*Pi);

a[k] := (1/Pi)*(Int(f(x)*cos(k*x), x = 0 .. 2*Pi))

 (10)
> 

a[0]:=subs(k=0,%):
> 

a[0]:=subs(cos(0)=1,%);

a[0] := (1/Pi)*(Int(f(x), x = 0 .. 2*Pi))

 (11)
> 

b[k]:=(1/Pi)*Int(f(x)*sin(k*x),x=0..2*Pi);

b[k] := (1/Pi)*(Int(f(x)*sin(k*x), x = 0 .. 2*Pi))

 (12)
> 

F(x):=Sp(x):
> 

A[0]:=value(subs(f(x)=F(x),a[0]));

A[0] := .9238343800

 (13)
> 

A[k]:=simplify(value(subs(f(x)=F(x),a[k]))):
> 

A[k]:=subs({sin(k*Pi)=0.,(cos(k*Pi))^2=1.},%):
> 

A[k]:=evalf(%):
> 

for i in [seq(k,k=1..5)] do A[i]:=subs(k=i,A[k]) od:
> 

for i in [seq(k,k=1..15)] do P[i]:=subs(k=i,A[k]) od:
> 

for i in [seq(k,k=1..100)] do R[i]:=subs(k=i,A[k]) od:
> 

B[k]:=simplify(value(subs(f(x)=F(x),b[k]))):
> 

B[k]:=subs({sin(k*Pi)=0.,(cos(k*Pi))^2=1.},%):
> 

for i in [seq(k,k=1..5)] do B[i]:=subs(k=i,B[k]) od:
> 

for i in [seq(k,k=1..15)] do Q[i]:=subs(k=i,B[k]) od:
> 

for i in [seq(k,k=1..100)] do S[i]:=subs(k=i,B[k]) od:
> 

# The coefficients A[i]...S[i] are not printed because of their length. However, it can be done by inserting semicolons (;) instead of doppelpoints (:) at the ends of the comands. The several FOURIER-Series with k = 5, 15, 100 are expressed, respectively, by Y(x), Z(x), T(x) as follows:

> 

 

> 

Y(x):=evalf(A[0]/2+sum(A[k]*cos(k*x)+B[k]*sin(k*x),k=1..5));

Y(x) := .4619171900-0.126839e-2*cos(x)-.1994095199*sin(x)-0.2112290000e-2*cos(2.*x)-0.6689752450e-1*sin(2.*x)+0.2246139506e-1*cos(3.*x)-0.7400420268e-1*sin(3.*x)+0.5560480578e-1*cos(4.*x)-0.3867736827e-1*sin(4.*x)-0.6746665347e-2*cos(5.*x)+0.3879916000e-2*sin(5.*x)

 (14)
> 

Z(x):=evalf(A[0]/2+sum(P[k]*cos(k*x)+Q[k]*sin(k*x),k=1..15));

Z(x) := 0.387991576e-2*sin(5.*x)-0.6746665105e-2*cos(5.*x)-0.3867736861e-1*sin(4.*x)+0.556048050e-1*cos(4.*x)-0.740042028e-1*sin(3.*x)+0.224613933e-1*cos(3.*x)-0.668975234e-1*sin(2.*x)-0.2112284e-2*cos(2.*x)-0.1268369784e-2*cos(x)-.1994095720*sin(x)+.4619171900+0.2949052345e-2*sin(14.*x)+0.1792042605e-1*cos(14.*x)-0.1146441471e-1*sin(13.*x)+0.6900361173e-2*cos(13.*x)+0.1562043141e-2*sin(12.*x)+0.9029515667e-2*cos(12.*x)-0.6726948932e-2*sin(11.*x)-0.8730614792e-2*cos(11.*x)+0.1098689589e-1*sin(10.*x)-0.1324355714e-1*sin(15.*x)-0.1926639980e-1*cos(15.*x)+0.3359667303e-1*cos(10.*x)-0.7154219204e-1*sin(9.*x)-0.7604331623e-2*cos(9.*x)-0.2814642451e-1*sin(8.*x)+0.1793538830e-1*cos(8.*x)-0.1755966185e-1*sin(7.*x)+0.6262061080e-2*cos(7.*x)-0.2997603416e-1*sin(6.*x)+0.1299643979e-1*cos(6.*x)

 (15)
> 

T(x):=evalf(A[0]/2+sum(R[k]*cos(k*x)+S[k]*sin(k*x),k=1..100)):
> 

# This equation is not printed because of its length with more than 100 terms. However, printing is possible,if we insert a semicolon (;) instead of the doppel point (:). The following Figure illustrades the approximations of the spline function by FOURIER-Series with  k = [5, 15, 100]):
> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35,                             title="Global Temperature Index (°C)"):

> 

p[3]:=plot({Y(x),Z(x),T(x)},x=0..2*Pi,th=2,co=black):
> 

p[4]:=plots[textplot]({[1.5,0.8,`DATA by Nasa`],                  [1.5,0.7,`Cubic Spline`],                                [4.5,0.25,`FOURIER-Series k = [5, 15, 100]`]}):
> 

plots[display](seq(p[k],k=1..4));

 

 

In this Figure we see that the FOURIER-Series with  k = [5,  15],  i.e., the functions  Y(x)  and  Z(x)  are not suitable approximations to the given DATA, because they do not increase in the period  x = [5.5.. 6.28] = [2010..2015]  like the cubic spline. The best approximation is the series with  k = 100, which is very close to the cubic spline. With increasing the value of  k one can improve the approximation. Finaly, with  k = infinity, we arrive exactly at the splinefunction.

 

> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35,                             title="Global Temperuture Index (°C)"):
> 

p[3]:=plot(T(x),x=0..2*Pi,th=2,co=black):
> 

p[4]:=plots[textplot]({[1.5,0.8,`DATA by Nasa`],                  [1.5,0.7,`Cubic Spline`],                                [4.5,0.25,`Fourier-Series  T(x)  with  k = 100`]}):
> 

plots[display](seq(p[k],k=1..4));

 

 

In this Figure we cannot observe a difference between the FOURIER-Series (k = 100) and the Cubic Spline.

The L[2] error norms can be expressed as:

 

> 

L[2][k=5]:=                               sqrt((1/Pi)*Int((Sp(xi)-Y(xi))^2,xi=Pi..2*Pi))= sqrt((1/Pi)*int((Sp(x)-Y(x))^2,x=Pi..2*Pi));

L[2][k = 5] := ((1/Pi)*(Int((Sp(xi)-Y(xi))^2, xi = Pi .. 2*Pi)))^(1/2) = 0.8181416648e-1

 (16)
> 

L[2][k=15]:=                      sqrt((1/Pi)*Int((Sp(xi)-Z(xi))^2,xi=Pi..2*Pi))= sqrt((1/Pi)*int((Sp(x)-Z(x))^2,x=Pi..2*Pi));

L[2][k = 15] := ((1/Pi)*(Int((Sp(xi)-Z(xi))^2, xi = Pi .. 2*Pi)))^(1/2) = 0.2918338151e-1

 (17)
> 

L[2][k=100]:=                     sqrt((1/Pi)*Int((Sp(xi)-T(xi))^2,xi=Pi..2*Pi))=0.01278068452;

L[2][k = 100] := ((1/Pi)*(Int((Sp(xi)-T(xi))^2, xi = Pi .. 2*Pi)))^(1/2) = 0.1278068452e-1

 (18)
> 

 

 

The three values above show the L[2] error norms of  the FOURIER-Series  Y(x),  Z(x),  and  T(x)  relative to the cubic spline function. Furthermore one should discuss the error norms with respect to the given experimental data, as follows:

 

> 

with(linalg):
> 

for i from 1 to 36 do                    u[i]:=subs(x=DATA[i][1],Y(x))-DATA[i][2] od:
> 

U:=evalf(vector([seq(u[i],i=1..36)]),4);

U := Vector[row](36, {(1) = .2498, (2) = 0.6659e-1, (3) = .1357, (4) = -.1408, (5) = -0.3049e-1, (6) = 0.3266e-1, (7) = 0.3307e-1, (8) = -0.3157e-1, (9) = -0.2774e-1, (10) = .1102, (11) = -0.5623e-1, (12) = -0.9228e-1, (13) = 0.6943e-1, (14) = 0.4757e-1, (15) = -0.4570e-2, (16) = -0.8647e-1, (17) = 0.8845e-1, (18) = 0.6887e-2, (19) = -.1322, (20) = 0.8777e-1, (21) = 0.8558e-1, (22) = -0.3127e-1, (23) = -0.7473e-1, (24) = -0.2643e-1, (25) = 0.8638e-1, (26) = -0.3500e-1, (27) = 0.1490e-1, (28) = -0.3613e-1, (29) = 0.6447e-1, (30) = -0.4337e-1, (31) = -0.9017e-1, (32) = 0.6089e-1, (33) = 0.7040e-1, (34) = 0.3150e-1, (35) = -.1087, (36) = -.3180})

 (19)
> 

l[2][k=5]:=                                 (1/sqrt(number_of_points))*Norm(U,2)=            evalf((1/sqrt(36))*norm(U,2),4);

l[2][k = 5] := Norm(U, 2)/number_of_points^(1/2) = 0.9762e-1

 (20)
> 

for i from 1 to 36 do                    v[i]:=subs(x=DATA[i][1],Z(x))-DATA[i][2] od:
> 

V:=evalf(vector([seq(v[i],i=1..36)]),4);

V := Vector[row](36, {(1) = .3189, (2) = -0.822e-1, (3) = 0.593e-1, (4) = -0.449e-1, (5) = 0.308e-1, (6) = -0.179e-1, (7) = 0.89e-2, (8) = -0.1e-3, (9) = -0.907e-2, (10) = 0.131e-1, (11) = -0.1048e-1, (12) = 0.222e-2, (13) = 0.1e-3, (14) = -0.17e-2, (15) = 0.51e-2, (16) = -0.8643e-2, (17) = 0.84e-2, (18) = -0.375e-2, (19) = 0.7e-3, (20) = -0.197e-2, (21) = 0.754e-2, (22) = -0.1251e-1, (23) = 0.163e-1, (24) = -0.136e-1, (25) = 0.1712e-1, (26) = -0.185e-1, (27) = 0.182e-1, (28) = -0.142e-1, (29) = 0.834e-2, (30) = 0.11e-2, (31) = -0.41e-2, (32) = 0.67e-2, (33) = -0.54e-2, (34) = -0.18e-2, (35) = 0.292e-1, (36) = -.2444})

 (21)
> 

l[2][k=15]:=                                (1/sqrt(number_of_points))*Norm(V,2)=            evalf((1/sqrt(36))*norm(V,2),4);

l[2][k = 15] := Norm(V, 2)/number_of_points^(1/2) = 0.7042e-1

 (22)
> 

for i from 1 to 36 do                    w[i]:=subs(x=DATA[i][1],T(x))-DATA[i][2] od:
> 

W:=evalf(vector([seq(w[i],i=1..36)]),4);

W := Vector[row](36, {(1) = .2869, (2) = 0.60e-2, (3) = -0.10e-2, (4) = -0.29e-2, (5) = -0.22e-2, (6) = -0.15e-2, (7) = 0., (8) = 0.11e-2, (9) = 0.127e-2, (10) = 0., (11) = -0.70e-3, (12) = -0.37e-3, (13) = -0.7e-3, (14) = -0.2e-3, (15) = 0.6e-3, (16) = 0.9498e-3, (17) = 0.12e-2, (18) = 0.19e-3, (19) = -0.4e-3, (20) = -0.36e-3, (21) = -0.16e-3, (22) = -0.70e-3, (23) = -0.3e-3, (24) = 0.9e-3, (25) = 0.17e-3, (26) = 0.7e-3, (27) = 0.13e-2, (28) = -0.6e-3, (29) = -0.194e-2, (30) = -0.6e-3, (31) = 0.3e-3, (32) = 0.27e-2, (33) = 0.28e-2, (34) = 0.22e-2, (35) = -0.89e-2, (36) = -.2270})

 (23)
> 

l[2][k=100]:=                               (1/sqrt(number_of_points))*Norm(W,2)=            evalf((1/sqrt(36))*norm(W,2),4);

l[2][k = 100] := Norm(W, 2)/number_of_points^(1/2) = 0.6102e-1

 (24)

 

 Note: Because of  its "piecevise character"  of  the spline function  the precision of  the L[2] error norm cannot be arrived by Digits, for instance, equal to four. Therefore, we have calculated the L[2] norms above with the default value of  ten. However, the calculation of  l[2] norms with Digits = 4  is possible without loss of  precision.

The error norms may sometimes be improved without changing the number of  terms, if  we utilize  the method of  smoothing. As an example let us apply this method to the FOURIER-Series  Y(x)  with  k = 5  terms. Introducing the smoothing factor  g(k,N) we arrive at the modified FOURIER-Siries  G(x,n,N) as:

 

> 

g(k,N):=N*sin(Pi*k/N)/Pi/k;   # smoothing factor

 

 

g(k, N) := N*sin(Pi*k/N)/(Pi*k)

 (25)
> 

G(x,n,N):= A[zero]/2+sum(g(kappa,Nu)*(A[kappa]*cos(kappa*x)+ B[kappa]*sin(kappa*x)),kappa=1..n);   # smoothing function

G(x, n, N) := (1/2)*A[zero]+sum(g(kappa, Nu)*(A[kappa]*cos(kappa*x)+B[kappa]*sin(kappa*x)), kappa = 1 .. n)

 (26)
> 

# Example for Y(x) with n = 5 and N = n + 1:
> 

g(k,6):=subs(N=6,g(k,N));

g(k, 6) := 6*sin(Pi*k/6)/(Pi*k)

 (27)
> 

G(x,n=5,N=6):=evalf(subs({n=5,Nu=6,kappa=k,g(kappa,Nu)=g(k,6), A[zero]=A[0],A[kappa]=A[k],B[kappa]=B[k]},G(x,n,N)));

G(x, n = 5, N = 6) := .4619171900-0.1211159545e-2*cos(x)-.1904220876*sin(x)-0.1746849536e-2*cos(2.*x)-0.5532380917e-1*sin(2.*x)+0.1429936731e-1*cos(3.*x)-0.4711253834e-1*sin(3.*x)+0.2299240155e-1*cos(4.*x)-0.1599296323e-1*sin(4.*x)-0.1288518137e-2*cos(5.*x)+0.7410093358e-3*sin(5.*x)

 (28)
> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35):
> 

p[3]:=plot({Y(x),G(x,n=5,N=6)},x=0..2*Pi,th=1,co=black,          title="Approximation with n = 5  and its Smoothing"):
> 

p[4]:=plots[textplot]({[2,0.75,`DATA by Nasa`],                    [2,0.65,`Cubic Spline`]}):
> 

plots[display](seq(p[k],k=1..4));

 
> 

# In this Figure both the FOURIER approximation Y(x) with n = 5 and its smoothing are represented by thin lines and differ from the cubic spline (thick line) in the period x = [5..2*Pi].The error Norm of the smoothing with respect to the experimental data is given as follows:

> 

for i from 1 to 36 do            m[i]:=subs(x=DATA[i][1],G(x,n=5,N=6))-DATA[i][2] od:
> 

M:=evalf(vector([seq(m[i],i=1..36)]),4);

M := Vector[row](36, {(1) = .2150, (2) = 0.7059e-1, (3) = .1803, (4) = -0.6808e-1, (5) = 0.4792e-1, (6) = 0.9122e-1, (7) = 0.5340e-1, (8) = -0.5277e-1, (9) = -0.7569e-1, (10) = 0.5742e-1, (11) = -0.8861e-1, (12) = -0.8804e-1, (13) = .1033, (14) = 0.9430e-1, (15) = 0.3175e-1, (16) = -0.7687e-1, (17) = 0.7087e-1, (18) = -0.2378e-1, (19) = -.1560, (20) = 0.8535e-1, (21) = .1008, (22) = -0.8413e-2, (23) = -0.6330e-1, (24) = -0.3211e-1, (25) = 0.6107e-1, (26) = -0.6475e-1, (27) = -0.7975e-3, (28) = -0.3190e-1, (29) = 0.8772e-1, (30) = -0.1613e-1, (31) = -0.7586e-1, (32) = 0.4567e-1, (33) = 0.2404e-1, (34) = -0.3168e-1, (35) = -.1693, (36) = -.3533})

 (29)
> 

l[2][G_k=5]:=                                 (1/sqrt(number_of_points))*Norm(M,2)=               evalf((1/sqrt(36))*norm(M,2),4);

l[2][G_k = 5] := Norm(M, 2)/number_of_points^(1/2) = .1038

 (30)

 

 

> 

# Another Example,i.e., the smoothing of the FOURIER-Series Z(x) with k = 15 and N = 16  has been illustrated, as follows:
> 

g(k,16):=subs(N=16,g(k,N));   # smoothing factor

g(k, 16) := 16*sin(Pi*k/16)/(Pi*k)

 (31)
> 

G(x,n=15,N=16):=               evalf(subs({n=15,Nu=16,kappa=k,g(kappa,Nu)=g(k,6), A[zero]=A[0],A[kappa]=A[k],B[kappa]=B[k]},G(x,n,N)));

G(x, n = 15, N = 16) := 0.7410093308e-3*sin(5.*x)-0.1288518136e-2*cos(5.*x)-0.1599296325e-1*sin(4.*x)+0.2299240155e-1*cos(4.*x)-0.4711253825e-1*sin(3.*x)+0.1429936729e-1*cos(3.*x)-0.5532380937e-1*sin(2.*x)-0.1746849128e-2*cos(2.*x)-0.1211174497e-2*cos(x)-.1904220903*sin(x)+.4619171900+0.3484066628e-3*sin(14.*x)+0.2117153272e-2*cos(14.*x)-0.8421315087e-3*sin(13.*x)+0.5068738156e-3*cos(13.*x)+0.5839784577e-3*sin(11.*x)+0.7579202581e-3*cos(11.*x)-0.1817217951e-2*sin(10.*x)-0.1686222065e-2*sin(15.*x)-0.2453074210e-2*cos(15.*x)-0.5556844983e-2*cos(10.*x)+0.1518172469e-1*sin(9.*x)+0.1613689340e-2*cos(9.*x)+0.5819226452e-2*sin(8.*x)-0.3708111709e-2*cos(8.*x)+0.2395463120e-2*sin(7.*x)-0.8542610962e-3*cos(7.*x)

 (32)
> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35):
> 

p[3]:=plot(Z(x),x=0..2*Pi,th=1,co=black,                    title="FOURIER-Series with k=15 and its Smoothing"):
> 

p[4]:=plot(G(x,n=15,N=16),x=0..2*Pi,th=2,co=black):
> 

p[5]:=plots[textplot]({[2,0.75,`DATA by Nasa`],                    [2,0.65,`Cubic Spline`]}):
> 

plots[display](seq(p[k],k=1..5));

 
> 

 

> 

# Error Norm of the smoothing with respect to the experimental data:
> 

for i from 1 to 36 do          q[i]:=subs(x=DATA[i][1],G(x,n=15,N=16))-DATA[i][2] od:
> 

Q:=evalf(vector([seq(q[i],i=1..36)]),4);

Q := Vector[row](36, {(1) = .2074, (2) = 0.9166e-1, (3) = .1904, (4) = -0.892e-1, (5) = 0.390e-1, (6) = .1107, (7) = 0.599e-1, (8) = -0.718e-1, (9) = -0.7558e-1, (10) = 0.695e-1, (11) = -0.8998e-1, (12) = -.1003, (13) = .1086, (14) = .1100, (15) = 0.129e-1, (16) = -0.8104e-1, (17) = 0.8700e-1, (18) = -0.2201e-1, (19) = -.1698, (20) = 0.8279e-1, (21) = .1213, (22) = -0.1499e-1, (23) = -0.7510e-1, (24) = -0.276e-1, (25) = 0.7182e-1, (26) = -0.671e-1, (27) = -0.170e-1, (28) = -0.253e-1, (29) = .1006, (30) = -0.172e-1, (31) = -0.907e-1, (32) = 0.397e-1, (33) = 0.478e-1, (34) = -0.305e-1, (35) = -.1896, (36) = -.3615})

 (33)
> 

l[2][G_k=15]:=                              (1/sqrt(number-of-points))*Norm(Q,2)=             evalf((1/sqrt(36))*norm(Q,2),4);

l[2][G_k = 15] := Norm(Q, 2)/(number-of-points)^(1/2) = .1110

 (34)

 

All Error norms calculated above proof that we have found suitable approximations to the experimental data. The following part of  this Worksheet is concerned with a simple model function, the parameters of  which have been determined by the

 

Nonlinear Regression with Maple

 

> 

restart:
> 

with(Statistics):   with(CurveFitting):
> 

# The experimental data are expressed as:
> 

x:=2*Pi*(xi-1980)/35;

x := 2*Pi*(xi-1980)/35

 (35)
> 

for i from 1980 to 2015 do x[i]:=[data]=2*Pi*(i-1980)/35 od:
> 

datax:=array([seq(2*Pi*(i-1980)/35,i=1980..2015)]);

datax := Vector[row](36, {(1) = 0, (2) = (2/35)*Pi, (3) = (4/35)*Pi, (4) = (6/35)*Pi, (5) = (8/35)*Pi, (6) = (2/7)*Pi, (7) = (12/35)*Pi, (8) = (2/5)*Pi, (9) = (16/35)*Pi, (10) = (18/35)*Pi, (11) = (4/7)*Pi, (12) = (22/35)*Pi, (13) = (24/35)*Pi, (14) = (26/35)*Pi, (15) = (4/5)*Pi, (16) = (6/7)*Pi, (17) = (32/35)*Pi, (18) = (34/35)*Pi, (19) = (36/35)*Pi, (20) = (38/35)*Pi, (21) = (8/7)*Pi, (22) = (6/5)*Pi, (23) = (44/35)*Pi, (24) = (46/35)*Pi, (25) = (48/35)*Pi, (26) = (10/7)*Pi, (27) = (52/35)*Pi, (28) = (54/35)*Pi, (29) = (8/5)*Pi, (30) = (58/35)*Pi, (31) = (12/7)*Pi, (32) = (62/35)*Pi, (33) = (64/35)*Pi, (34) = (66/35)*Pi, (35) = (68/35)*Pi, (36) = 2*Pi})

 (36)
> 

X:=evalf(array([seq(2*Pi*(i-1980)/35,i=1980..2015)]),3);

X := Vector[row](36, {(1) = 0., (2) = .179, (3) = .358, (4) = .537, (5) = .719, (6) = .898, (7) = 1.08, (8) = 1.26, (9) = 1.43, (10) = 1.61, (11) = 1.79, (12) = 1.98, (13) = 2.15, (14) = 2.33, (15) = 2.51, (16) = 2.69, (17) = 2.87, (18) = 3.05, (19) = 3.23, (20) = 3.42, (21) = 3.58, (22) = 3.77, (23) = 3.96, (24) = 4.11, (25) = 4.30, (26) = 4.49, (27) = 4.68, (28) = 4.84, (29) = 5.02, (30) = 5.21, (31) = 5.37, (32) = 5.56, (33) = 5.75, (34) = 5.93, (35) = 6.09, (36) = 6.28})

 (37)
> 

whattype(X);

symbol

 (38)
> 

Y:=array([0.28,0.33,0.13,0.31,0.16,0.12,0.19,0.34, 0.40,0.29,0.44,0.43,0.23,0.24,0.32,0.46,0.35,0.48,0.64, 0.42,0.42,0.55,0.63,0.62,0.55,0.69,0.63,0.66,0.54,0.65, 0.72,0.61,0.63,0.66,0.75,0.85]);

Y := Vector[row](36, {(1) = .28, (2) = .33, (3) = .13, (4) = .31, (5) = .16, (6) = .12, (7) = .19, (8) = .34, (9) = .40, (10) = .29, (11) = .44, (12) = .43, (13) = .23, (14) = .24, (15) = .32, (16) = .46, (17) = .35, (18) = .48, (19) = .64, (20) = .42, (21) = .42, (22) = .55, (23) = .63, (24) = .62, (25) = .55, (26) = .69, (27) = .63, (28) = .66, (29) = .54, (30) = .65, (31) = .72, (32) = .61, (33) = .63, (34) = .66, (35) = .75, (36) = .85})

 (39)
> 

whattype(Y);

symbol

 (40)
> 

DATA:=seq([X[i],Y[i]],i=1..36);

DATA := [0., .28], [.179, .33], [.358, .13], [.537, .31], [.719, .16], [.898, .12], [1.08, .19], [1.26, .34], [1.43, .40], [1.61, .29], [1.79, .44], [1.98, .43], [2.15, .23], [2.33, .24], [2.51, .32], [2.69, .46], [2.87, .35], [3.05, .48], [3.23, .64], [3.42, .42], [3.58, .42], [3.77, .55], [3.96, .63], [4.11, .62], [4.30, .55], [4.49, .69], [4.68, .63], [4.84, .66], [5.02, .54], [5.21, .65], [5.37, .72], [5.56, .61], [5.75, .63], [5.93, .66], [6.09, .75], [6.28, .85]

 (41)
> 

# The cubic splinefunction is given by:
> 

Sp(x):=Spline([DATA],x,degree=3):
> 

# A suitable model function to the DATA has been assumed as:
> 

f:=(t,A,B,a,b,c,d,e) -> a+A*exp(B*t)*(b*cos(c*t)+d*sin(e*t));

f := proc (t, A, B, a, b, c, d, e) options operator, arrow; a+A*exp(B*t)*(b*cos(c*t)+d*sin(e*t)) end proc

 (42)
> 

y:=unapply(evalf(NonlinearFit(f(t,A,B,a,b,c,d,e),X,Y,t),4),t);

y := proc (t) options operator, arrow; .3148+0.3582e-2*exp(.9353*t)*(-2.897*cos(1.239*t)-3.728*sin(.9736*t)) end proc

 (43)
> 

y(x):=subs(t=x,y(t));

y(x) := .3148+0.3582e-2*exp(.9353*x)*(-2.897*cos(1.239*x)-3.728*sin(.9736*x))

 (44)
> 

y(0):=simplify(subs(x=0,y(x)));

y(0) := .304422946

 (45)
> 

y(0)[data]:=0.28;

y(0)[data] := .28

 (46)
> 

y(2*Pi):=simplify(subs(x=2*Pi,y(x)));

y(2*Pi) := .8455447604

 (47)
> 

y(2*Pi)[data]:=0.85;

y(2*Pi)[data] := .85

 (48)
> 

# The results are representd in the next Figure:
> 

 

> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35,                             title="Global Temperature Index (°C)"):
> 

p[3]:=plot(y(x),x=0..2*Pi,th=3,co=black):
> 

p[4]:=plots[textplot]({[2,0.75,`DATA by Nasa`],                    [2,0.65,`Cubic Spline`],[4.5,0.3,`Nonlinear Regression`], [4.5,0.2,`Approximation y(x)`]}):
> 

plots[display](seq(p[k],k=1..4));

 

 

In this Figure, we see, the approximation  y(x) is verry close to the cubic spline.  The L[2] error norm between the cubic spline and the approximation  y(x)  can be expressed as:

 

> 

L[2]:=                         sqrt((1/(2*Pi))*Int((SPLINE-APPROXIMATION)^2,xi=0..2*Pi))= sqrt((1/2/Pi)*int((Sp(x)-y(x))^2,x=0..2*Pi));

L[2] := (1/2)*(2^(1/2)*((1/Pi)*(Int((SPLINE-APPROXIMATION)^2, xi = 0 .. 2*Pi)))^(1/2)) = 0.8181141513e-1

 (49)

The l[2] error norm of  the approximation  y(x)  with respect to the experimental data is given by:

 

> 

with(linalg):
> 

for i from 1 to 36 do                    u[i]:=subs(x=DATA[i][1],y(x))-DATA[i][2] od:
> 

U:=vector([seq(u[i],i=1..36)]);

U := Vector[row](36, {(1) = 0.24422946e-1, (2) = -0.2990521279e-1, (3) = .1653248282, (4) = -0.1970833595e-1, (5) = .1251637815, (6) = .1604439126, (7) = 0.8639540586e-1, (8) = -0.6637435981e-1, (9) = -.1273699943, (10) = -0.1615192083e-1, (11) = -.1620153257, (12) = -.1438610605, (13) = 0.6719203757e-1, (14) = 0.7311090985e-1, (15) = 0.1353535615e-1, (16) = -.1015990652, (17) = 0.3735430600e-1, (18) = -0.602825065e-1, (19) = -.1855268186, (20) = 0.722679495e-1, (21) = .1036078669, (22) = 0.84059714e-2, (23) = -0.414537161e-1, (24) = -0.122427551e-1, (25) = 0.750205663e-1, (26) = -0.563249393e-1, (27) = 0.40344509e-2, (28) = -0.308672925e-1, (29) = 0.804514501e-1, (30) = -0.380588973e-1, (31) = -.1103774758, (32) = 0.85249307e-2, (33) = 0.162382905e-1, (34) = 0.340846772e-1, (35) = 0.47840012e-2, (36) = -0.60401312e-2})

 (50)
> 

l[2]:=                                      (1/sqrt(number_of_points))*Norm(U,2)=                      evalf((1/sqrt(36))*norm(U,2),4);

l[2] := Norm(U, 2)/number_of_points^(1/2) = 0.8415e-1

 (51)

 

Both the L[2] and the l[2] error norms show that the assumed function  y(x)  is a suitable approximation to the experimental data.

Note:  Because of  its  "piecevise character"  of  the splinefunction the high precision of  the  L[2] norm cannot be arrived by Digits, for instance, equal to four. Thus, we have calculated the L[2] error norms above with "default" value of ten. However, the calculation of  l[2] error norms with Digits = 4 is possible without loss of  precision.

   Another example within the category  Nonlinear Regression  should be discussed in the following:

 

> 

h:=(t,A,B,C,D,a,b,c) -> a+A*exp(B*t)*cos(b*t)+C*exp(D*t)*sin(c*t);

 

 

h := proc (t, A, B, C, D, a, b, c) options operator, arrow; a+A*exp(B*t)*cos(b*t)+C*exp(D*t)*sin(c*t) end proc

 (52)
> 

z:=unapply(evalf(NonlinearFit(h(t,A,B,C,D,a,b,c),X,Y,t),4),t);

z := proc (t) options operator, arrow; .2925-0.8539e-2*exp(.9705*t)*cos(1.001*t)-0.1584e-2*exp(1.267*t)*sin(.7948*t) end proc

 (53)
> 

z(x):=subs(t=x,z(t));

z(x) := .2925-0.8539e-2*exp(.9705*x)*cos(1.001*x)-0.1584e-2*exp(1.267*x)*sin(.7948*x)

 (54)
> 

z(0):=simplify(subs(x=0,z(x)));

z(0) := .283961

 (55)
> 

z(0)[data]:=0.28;

z(0)[data] := .28

 (56)
> 

z(2*Pi):=simplify(subs(x=2*Pi,z(x)));

z(2*Pi) := .855255878

 (57)
> 

z(2*Pi)[data]:=0.85;

z(2*Pi)[data] := .85

 (58)
> 

# The results of the approximation z(x) are represented in the next Figure:
> 

alias(th=thickness,co=color):
> 

p[1]:=plot(Sp(x),x=0..2*Pi,axes=boxed,th=3,co=black):
> 

p[2]:=plot([DATA],x=0..2*Pi,th=3,co=black, style=point,symbol=cross,symbolsize=35,                             title="Global Temperature Index (°C)"):
> 

p[3]:=plot(z(x),x=0..2*Pi,th=3,co=black):
> 

p[4]:=plots[textplot]({[2,0.75,`DATA by Nasa`],           [2,0.65,`CubicSpline`], [4.5,0.3,`Nonlinear Regression`], [4.5,0.2,`Approximation z(x)`]}):
> 

plots[display](seq(p[k],k=1..4));

 

 

In this Figure, we see, the approximation  z(x)  is very close to the cubic spline, similar to the approimation  y(x)  in one Figure before. The M[2] error norm between the cubic spline and the approximation  z(x)  can be expressed as:

 

> 

M[2]:=                   sqrt((1/2/Pi)*Int((SPLINE-APPROXIMATION)^2,xi=0..2*Pi))= sqrt((1/2/Pi)*int((Sp(x)-z(x))^2,x=0..2*Pi));

M[2] := (1/2)*(2^(1/2)*((1/Pi)*(Int((SPLINE-APPROXIMATION)^2, xi = 0 .. 2*Pi)))^(1/2)) = 0.8037261461e-1

 (59)

 

The m[2] error norm of  the approximation  z(x)  with respect to the experimatal data is given by:

 

> 

with(linalg):
> 

for i from 1 to 36 do                    v[i]:=subs(x=DATA[i][1],z(x))-DATA[i][2] od:
> 

V:=vector([seq(v[i],i=1..36)]);

V := Vector[row](36, {(1) = 0.3961e-2, (2) = -0.4777816197e-1, (3) = .1504815197, (4) = -0.3114638782e-1, (5) = .1174672875, (6) = .1565585479, (7) = 0.8633369230e-1, (8) = -0.6291983162e-1, (9) = -.1210486070, (10) = -0.7505475104e-2, (11) = -.1519988774, (12) = -.1336479029, (13) = 0.7636888008e-1, (14) = 0.7989929510e-1, (15) = 0.1667348684e-1, (16) = -.1031763805, (17) = 0.3033314104e-1, (18) = -0.7300329763e-1, (19) = -.2036168606, (20) = 0.4959397336e-1, (21) = 0.7849631705e-1, (22) = -0.1732185911e-1, (23) = -0.6490280976e-1, (24) = -0.3178133171e-1, (25) = 0.627937730e-1, (26) = -0.594905968e-1, (27) = 0.101314769e-1, (28) = -0.183286679e-1, (29) = 0.969226508e-1, (30) = -0.231649540e-1, (31) = -.102295258, (32) = 0.2491299e-2, (33) = -0.6965774e-2, (34) = -0.732068e-3, (35) = -0.26510733e-1, (36) = 0.2523863e-2})

 (60)
> 

m[2]:=                                      (1/sqrt(number_of_points))*Norm(V,2)=             evalf((1/sqrt(36))*norm(V,2),4);

m[2] := Norm(V, 2)/number_of_points^(1/2) = 0.8255e-1

 (61)

 

Both the M[2]  and  the m[2]  error norms show that the assumed function  z(x) on pages 14/15 like  y(x) on pages 12/13 is a suitable approximation to the experimental data. The error norms   M[2], m[2]   of  the last approximation  z(x)  are smaler, respectively, than the error norms   L[2], l[2]   of  the approximation  y(x).

Note:   Because of its  "piecevise character"  of  the splinefunction the high precision of  the M[2] error norm cannot be arrived by Digits, for example, equal to four. Thus, we have calculated the M[2] error norm with  "default"  value of  ten. However, the calculation of  the m[2] error norm with  Digits = 4  is possible without loss of  precision.

 

Five-Year Mean Temperature Change with the Base Period  1950...2015

 

> 

restart:
> 

with(Statistics):  with(CurveFitting):
> 

dataxl:=[seq(k,k=1950..2015,5)];

dataxl := [1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015]

 (62)
> 

whattype(dataxl);

list

 (63)
> 

X:=array([seq(k,k=1950..2015,5)]);

X := Vector[row](14, {(1) = 1950, (2) = 1955, (3) = 1960, (4) = 1965, (5) = 1970, (6) = 1975, (7) = 1980, (8) = 1985, (9) = 1990, (10) = 1995, (11) = 2000, (12) = 2005, (13) = 2010, (14) = 2015})

 (64)
> 

whattype(X);

symbol

 (65)
> 

datay[1946..1950]:=[-0.04,-0.04,-0.10,-0.10,-0.18];

datay[1946 .. 1950] := [-0.4e-1, -0.4e-1, -.10, -.10, -.18]

 (66)
> 

m[1950]:=evalf(Mean(%),3);

m[1950] := -0.920e-1

 (67)
> 

datay[1951..1955]:=[-0.07,0.01,0.09,-0.12,-0.14];

datay[1951 .. 1955] := [-0.7e-1, 0.1e-1, 0.9e-1, -.12, -.14]

 (68)
> 

m[1955]:=evalf(Mean(%),3);

m[1955] := -0.460e-1

 (69)
> 

datay[1956..1960]:=[-0.19,0.04,0.06,0.03,-0.03];

datay[1956 .. 1960] := [-.19, 0.4e-1, 0.6e-1, 0.3e-1, -0.3e-1]

 (70)
> 

m[1960]:=evalf(Mean(%),3);

m[1960] := -0.180e-1

 (71)
> 

datay[1961..1965]:=[0.05,0.03,0.06,-0.21,-0.10];

datay[1961 .. 1965] := [0.5e-1, 0.3e-1, 0.6e-1, -.21, -.10]

 (72)
> 

m[1965]:=evalf(Mean(%),3);

m[1965] := -0.340e-1

 (73)
> 

datay[1966..1970]:=[-0.05,-0.02,-0.07,0.06,0.03];

datay[1966 .. 1970] := [-0.5e-1, -0.2e-1, -0.7e-1, 0.6e-1, 0.3e-1]

 (74)
> 

m[1970]:=evalf(Mean(%),3);

m[1970] := -0.100e-1

 (75)
> 

data[1972..1975]:=[-0.09,0.01,0.15,-0.08,-0.01];

data[1972 .. 1975] := [-0.9e-1, 0.1e-1, .15, -0.8e-1, -0.1e-1]

 (76)
> 

m[1975]:=evalf(Mean(%),3);

m[1975] := -0.400e-2

 (77)
> 

datay[1976..1980]:=[-0.11,0.18,0.07,0.17,0.28];

datay[1976 .. 1980] := [-.11, .18, 0.7e-1, .17, .28]

 (78)
> 

m[1980]:=evalf(Mean(%),3);

m[1980] := .118

 (79)
> 

datay[1981..1985]:=[0.33,0.13,0.31,0.16,0.12];

datay[1981 .. 1985] := [.33, .13, .31, .16, .12]

 (80)
> 

m[1985]:=evalf(Mean(%),3);

m[1985] := .210

 (81)
> 

datay[1986..1990]:=[0.19,0.34,0.40,0.29,0.44];

datay[1986 .. 1990] := [.19, .34, .40, .29, .44]

 (82)
> 

m[1990]:=evalf(Mean(%),4);

m[1990] := .3320

 (83)
> 

datay[1991..1995]:=[0.43,0.23,0.24,0.32,0.46];

datay[1991 .. 1995] := [.43, .23, .24, .32, .46]

 (84)
> 

m[1995]:=evalf(Mean(%),4);

m[1995] := .3360

 (85)
> 

datay[1996..2000]:=[0.35,0.48,0.64,0.42,0.42];

datay[1996 .. 2000] := [.35, .48, .64, .42, .42]

 (86)
> 

m[2000]:=evalf(Mean(%),4);

m[2000] := .4620

 (87)
> 

datay[2001..2005]:=[0.55,0.63,0.62,0.55,0.69];

datay[2001 .. 2005] := [.55, .63, .62, .55, .69]

 (88)
> 

m[2005]:=evalf(Mean(%),4);

m[2005] := .6080

 (89)
> 

datay[2006..2010]:=[0.63,0.66,0.54,0.65,0.72];

datay[2006 .. 2010] := [.63, .66, .54, .65, .72]

 (90)
> 

m[2010]:=evalf(Mean(%),3);

m[2010] := .640

 (91)
> 

datay[2011..2015]:=[0.61,0.63,0.66,0.75,0.85];

datay[2011 .. 2015] := [.61, .63, .66, .75, .85]

 (92)
> 

m[2015]:=evalf(Mean(%),3);

m[2015] := .700

 (93)
> 

datayl:=[seq(m[i],i=1950..2015,5)];

datayl := [-0.920e-1, -0.460e-1, -0.180e-1, -0.340e-1, -0.100e-1, -0.400e-2, .118, .210, .3320, .3360, .4620, .6080, .640, .700]

 (94)
> 

whattype(datayl);

list

 (95)
> 

Y:=array([seq(m[i],i=1950..2015,5)]);

Y := Vector[row](14, {(1) = -0.920e-1, (2) = -0.460e-1, (3) = -0.180e-1, (4) = -0.340e-1, (5) = -0.100e-1, (6) = -0.400e-2, (7) = .118, (8) = .210, (9) = .3320, (10) = .3360, (11) = .4620, (12) = .6080, (13) = .640, (14) = .700})

 (96)
> 

whattype(Y);

symbol

 (97)
> 

DATAL:=[seq([i,m[i]],i=1950..2015,5)];

DATAL := [[1950, -0.920e-1], [1955, -0.460e-1], [1960, -0.180e-1], [1965, -0.340e-1], [1970, -0.100e-1], [1975, -0.400e-2], [1980, .118], [1985, .210], [1990, .3320], [1995, .3360], [2000, .4620], [2005, .6080], [2010, .640], [2015, .700]]

 (98)
> 

whattype(DATAL);

list

 (99)
> 

DATA:=seq([i,m[i]],i=1950..2015,5);

DATA := [1950, -0.920e-1], [1955, -0.460e-1], [1960, -0.180e-1], [1965, -0.340e-1], [1970, -0.100e-1], [1975, -0.400e-2], [1980, .118], [1985, .210], [1990, .3320], [1995, .3360], [2000, .4620], [2005, .6080], [2010, .640], [2015, .700]

 (100)
> 

whattype(DATA);

exprseq

 (101)

Cubic Spline:

 

> 

Sp(x):=Spline([DATA],x,degree=3):
> 

alias(th=thickness,co=color):
> 

p[1]:=plot([DATA],x=1950..2015,axes=boxed,th=3,co=black, style=point,symbol=cross,symbolsize=30):
> 

p[2]:=plot(Sp(x),x=1950..2015,th=3,co=black,                     title="Five-Year Mean Temperature Index (°C)  Period 1950...2015"):
> 

p[3]:=plot(0,x=1950..2015,th=1,co=black):
> 

p[4]:=plots[textplot]({[1970,0.6,`DATA by NASA`],                [1970,0.4,`Cubic Spline Interpolation`]}):
> 

plots[display](seq(p[k],k=1..4));

 
> 

 

Nonlinear Regression:  Cubic Polynom

 

> 

with(Statistics):
> 

f:=(t,a,b)->-0.092+a*((t-1950)/65)+b*((t-1950)/65)^2+ (0.792-a-b)*((t-1950)/65)^3;

f := proc (t, a, b) options operator, arrow; -0.92e-1+a*((1/65)*t-30)+b*((1/65)*t-30)^2+(.792-a-b)*((1/65)*t-30)^3 end proc

 (102)
> 

y:=unapply(evalf(NonlinearFit(f(t,a,b),X,Y,t)),t);

y := proc (t) options operator, arrow; HFloat(5.228359901129614)-HFloat(0.0027283896928869816)*t+HFloat(1.8408343419295945)*(0.1538461538e-1*t-30.)^2-HFloat(0.8714890118919407)*(0.1538461538e-1*t-30.)^3 end proc

 (103)
> 

y(x):=subs(t=x,y(t));

y(x) := 5.22835990112961-0.272838969288698e-2*x+1.84083434192959*(0.1538461538e-1*x-30.)^2-.871489011891941*(0.1538461538e-1*x-30.)^3

 (104)
> 

y(1950):=evalf(subs(x=1950,y(x)));   y(1950)[data]:=-0.0920;

y(1950) := -0.920000000000003e-1

y(1950)[data] := -0.920e-1

 (105)
> 

y(2015):=evalf(subs(x=2015,y(x)));   y(2015)[data]:=0.700;

y(2015) := .699999989327983

y(2015)[data] := .700

 (106)

  The Boundary Conditions are fulfilled exactly.

 

> 

alias(th=thickness,co=color):
> 

p[1]:=plot([DATA],x=1950..2015,axes=boxed,th=3,co=black, style=point,symbol=cross,symbolsize=30):
> 

p[2]:=plot(Sp(x),x=1950..2015,th=3,co=black,                     title="Five-Year Mean Temperature Index (°C)  Period 1950...2015"):
> 

p[3]:=plot(0,x=1950..2015,th=1,co=black):
> 

p[4]:=plot(y(x),x=1950..2015,th=2,co=black):
> 

p[5]:=plots[textplot]({[1970,0.6,`DATA by NASA`],                 [1970,0.5,`Cubic Spline`],                                   [1970,0.4,`Nonlinear Regression`],                       [1970,0.3,`Approximation y(x)`]}):
> 

plots[display](seq(p[k],k=1..5));

 
> 

 

The L[2] error norm between the cubic spline and the approximation  y(x)  can be expressed as:

 

> 

L[2]:=                            sqrt((1/65)*Int((SPLINE-Z(x))^2,Z=1950..2015))= evalf(sqrt((1/65)*int((Sp(x)-y(x))^2,x=1950..2015)),4);

L[2] := (1/65)*(65^(1/2)*(Int((SPLINE-Z(x))^2, Z = 1950 .. 2015))^(1/2)) = 0.3621e-1

 (107)

 

The  l[2]  error norm of  the approximation  y(x)  with respect to the experimental data is given by:

 

> 

with(linalg):
> 

for i in [seq(1950..2015,5)] do v[i]:=subs(x=i,y(x))-m[i] od:
> 

V:=evalf(vector([seq(v[i],i=1950..2015,5)]),4);

V := Vector[row](14, {(1) = -0.3469000000e-15, (2) = -0.4915e-1, (3) = -0.6089e-1, (4) = -0.1160e-1, (5) = 0.1233e-1, (6) = 0.6652e-1, (7) = 0.1460e-1, (8) = 0.1808e-3, (9) = -0.3911e-1, (10) = 0.4234e-1, (11) = 0.2159e-2, (12) = -0.6004e-1, (13) = -0.1263e-1, (14) = -0.1067000000e-7})

 (108)
> 

l[2]:=                                      (1/sqrt(number_of_points))*Norm(V,2)=               evalf((1/sqrt(14))*norm(V,2),4);

l[2] := Norm(V, 2)/number_of_points^(1/2) = 0.3600e-1

 (109)
> 

 

> 

# Example ->  exp(...) and tanh(...)  combined:  
> 

h:=                                                                     (t,a,b)-> -0.092+   (0.792/(tanh(1)*exp(b)))*exp(b*(t-1950)/65)*tanh((t-1950)/65);

h := proc (t, a, b) options operator, arrow; -0.92e-1+.792*exp((1/65)*(b*(t-1950)))*tanh((1/65)*t-30)/((tanh(1)*exp(b))) end proc

 (110)
> 

r:=unapply(evalf(NonlinearFit(h(t,a,b),X,Y,t)),t);

r := proc (t) options operator, arrow; -0.92e-1+HFloat(0.3289842963030672)*exp(HFloat(0.017706043728021637)*t-HFloat(34.5267852696422))*tanh(0.1538461538e-1*t-30.) end proc

 (111)
> 

r(x):=subs(t=x,r(t));

r(x) := -0.92e-1+.328984296303067*exp(0.177060437280216e-1*x-34.5267852696422)*tanh(0.1538461538e-1*x-30.)

 (112)
> 

r(1950):=evalf(subs(x=1950,r(x)));  r(1950)[data]:=-0.0920;

r(1950) := -0.920000032898430e-1

r(1950)[data] := -0.920e-1

 (113)
> 

r(2015):=evalf(subs(x=2015,r(x)));  r(2015)[data]:=0.700;

r(2015) := .699999995279239

r(2015)[data] := .700

 (114)
> 

# The boundary conditions are exactly fulfilled.
> 

alias(th=thickness,co=color):
> 

p[1]:=plot([DATA],x=1050..2015,axes=boxed,th=3,co=black, style=point,symbol=cross,symbolsize=30):
> 

p[2]:=plot(Sp(x),x=1950..2015,th=3,co=black,                      title="Five-Year Mean Temperature Index (°C)  Period  1950...2015"):
> 

p[3]:=plot(0,x=1950...2015,th=1,co=black):
> 

p[4]:=plot(r(x),x=1950..2015,th=2,co=black):
> 

p[5]:=plots[textplot]({[1970,0.6,`DATA by NASA`],                [1970,0.5,`Cubic Spline`],                                   [1970,0.4,`Nonlinear Regression`],                       [1970,0.3,`Approximation  r(x)`]}):
> 

plots[display](seq(p[k],k=1..5));

 
> 

 

The  N[2]  error norm  between the cubic spline and the approximation  r(x)  can be expressed as:

 

> 

N[2]:=sqrt((1/65)*Int((SPLINE-R(x))^2,x=1950..2015))= evalf(sqrt((1/65)*int((Sp(x)-r(x))^2,x=1950..2015)));

N[2] := (1/65)*(65^(1/2)*(Int((SPLINE-R(x))^2, x = 1950 .. 2015))^(1/2)) = 0.4940573100e-1

 (115)
> 

for i in [seq(1950..2015,5)] do s[i]:=subs(x=i,r(x))-m[i] od:
> 

S:=evalf(vector([seq(s[i],i=1950..2015,5)]));

S := Vector[row](14, {(1) = -0.3289842963e-8, (2) = -0.184053718475058e-1, (3) = -0.140553325780236e-1, (4) = 0.392932303190068e-1, (5) = 0.578536888142986e-1, (6) = 0.998181630225897e-1, (7) = 0.313688694424441e-1, (8) = -0.130869094603342e-2, (9) = -0.580118795428445e-1, (10) = 0.949270969804600e-2, (11) = -0.385181108309168e-1, (12) = -0.997129538157974e-1, (13) = -0.396962268990830e-1, (14) = -0.4622868977e-8})

 (116)
> 

n[2]:=                                      (1/sqrt(number_of_points))*Norm(S,2)=               evalf((1/sqrt(14))*norm(S,2),4);

n[2] := Norm(S, 2)/number_of_points^(1/2) = 0.4842e-1

 (117)

 

 

Résumé:  The  error norms   relativ to the spline interpolations or alternatively with respect to  experimental data have

shown that the selected nonlinear model functions are suitable to fit the given data.  

 

>