tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

@Jsevillamol 

The attached works.

My original response was developed in Maple 2018

Note however that in Maple 2017 the *appears* to be sensitive to the location of the with(Domains) command - because the second execution group in the attached does not work!! Although it does in 2018 - see my earlier post. No idea why!!

(I fixed the GI/ED stuff as well)

restart;

with(Domains):
EEA:=proc(ED,a,b)
    description
    "Extended Euclidean Algorithm"
    "INPUT: an Euclidean Domain ED and two elements from said domain"
    "OUTPUT: r,s,t such that r = gcd(a,b) = s*a + t*b ";
    local r_0, r_1, r_aux, s_0, s_1, s_aux, t_0, t_1, t_aux, q;
    # Domain checks
    # TODO: check that ED is an euclidean domain
    if not ED[Type](a) then error "1st argument must be of type ED" end if;
    if not ED[Type](b) then error "2nd argument must be of type ED" end if;
    
    # Initialization
    r_0 := a; r_1 := b; # gcd series
    s_0 := ED[Input](1); s_1 := ED[Input](0); # 1st cofactor series
    t_0 := ED[Input](0); t_1 := ED[Input](1); # 2nd cofactor series

    # Loop

    while ED[`<>`](r_1,ED[Input](0)) do;
        print("All is fine before the Quo");
        print(r_0); print(r_1);
        q := ED[Quo](r_0, r_1);
        print("All is fine after the Quo");       
     #   r_aux := r_0 - q * r_1;
         r_aux:= ED[`-`](r_0, ED[`*`](q, r_1));
         r_0 := r_1; r_1 := r_aux;       
     #   s_aux := s_0 - q * s_1;
         s_aux:= ED[`-`](s_0, ED[`*`](q, s_1));
         s_0 := s_1; s_1 := s_aux;

     #   t_aux := t_0 - q * t_1;
         t_aux:= ED[`-`](t_0, ED[`*`](q, t_1));
         t_0 := t_1; t_1 := t_aux;
    od;
    # Result
    return r_0, s_0, t_0;
    end proc:
GI := Gaussian(Z):
a := GI[Input](-87+47*I);
b := GI[Input](-90+43*I);
r, s, t := EEA(GI, a, b);

---------------------- Domains version 1.0 ---------------------
Initially defined domains are Z and Q the integers and rationals

Abbreviations, e.g. DUP for DenseUnivariatePolynomial, also made

 

`domains/Gaussian/badge0`(-87, 47)

 

`domains/Gaussian/badge0`(-90, 43)

 

"All is fine before the Quo"

 

`domains/Gaussian/badge0`(-87, 47)

 

`domains/Gaussian/badge0`(-90, 43)

 

"All is fine after the Quo"

 

"All is fine before the Quo"

 

`domains/Gaussian/badge0`(-90, 43)

 

`domains/Gaussian/badge0`(3, 4)

 

"All is fine after the Quo"

 

"All is fine before the Quo"

 

`domains/Gaussian/badge0`(3, 4)

 

`domains/Gaussian/badge0`(2, -1)

 

"All is fine after the Quo"

 

"All is fine before the Quo"

 

`domains/Gaussian/badge0`(2, -1)

 

`domains/Gaussian/badge0`(1, 0)

 

"All is fine after the Quo"

 

`domains/Gaussian/badge0`(1, 0), `domains/Gaussian/badge0`(-39, -8), `domains/Gaussian/badge0`(39, 6)

(1)

restart;

EEA:=proc(ED,a,b)
    description
    "Extended Euclidean Algorithm"
    "INPUT: an Euclidean Domain ED and two elements from said domain"
    "OUTPUT: r,s,t such that r = gcd(a,b) = s*a + t*b ";
    local r_0, r_1, r_aux, s_0, s_1, s_aux, t_0, t_1, t_aux, q;
    # Domain checks
    # TODO: check that ED is an euclidean domain
    if not ED[Type](a) then error "1st argument must be of type ED" end if;
    if not ED[Type](b) then error "2nd argument must be of type ED" end if;
    
    # Initialization
    r_0 := a; r_1 := b; # gcd series
    s_0 := ED[Input](1); s_1 := ED[Input](0); # 1st cofactor series
    t_0 := ED[Input](0); t_1 := ED[Input](1); # 2nd cofactor series

    # Loop

    while ED[`<>`](r_1,ED[Input](0)) do;
        print("All is fine before the Quo");
        print(r_0); print(r_1);
        q := ED[Quo](r_0, r_1);
        print("All is fine after the Quo");       
     #   r_aux := r_0 - q * r_1;
         r_aux:= ED[`-`](r_0, ED[`*`](q, r_1));
         r_0 := r_1; r_1 := r_aux;       
     #   s_aux := s_0 - q * s_1;
         s_aux:= ED[`-`](s_0, ED[`*`](q, s_1));
         s_0 := s_1; s_1 := s_aux;

     #   t_aux := t_0 - q * t_1;
         t_aux:= ED[`-`](t_0, ED[`*`](q, t_1));
         t_0 := t_1; t_1 := t_aux;
    od;
    # Result
    return r_0, s_0, t_0;
    end proc:
with(Domains):
GI := Gaussian(Z):
a := GI[Input](-87+47*I);
b := GI[Input](-90+43*I);
r, s, t := EEA(GI, a, b);

---------------------- Domains version 1.0 ---------------------
Initially defined domains are Z and Q the integers and rationals

Abbreviations, e.g. DUP for DenseUnivariatePolynomial, also made

 

`domains/Gaussian/badge0`(-87, 47)

 

`domains/Gaussian/badge0`(-90, 43)

 

Error, (in EEA) cannot determine if this expression is true or false: not GI[Type](`domains/Gaussian/badge0`(-87, 47))

 

 


 

Download doms2017.mw

I assume that you are using the LinearAlgebra[Generic] package and have some modules somewhere which define various Euclidean domains. Without these, no-one can rerun/diagnose your problem: the only option is simple 'reading' of your code - which is a surprisingly difficult way to diagnose errors!

One thing which bothers me is that the two cases you cite are actually different - relating to representation of  the square root of -1. In the case you describe as "working", ie

GI:=Gaussians(Z): a := GI[Input](-90+43*I); b := GI[Input](3+4*I); GI[Quo](a, b)

you have the conventional 'I' in the input, although '_i' appears in the output.

On the other hand in the procedure, you always seem to have _i - both input and output?? Presumably when you first call this procedure, both the arguments 'a' and 'b' actuall have the the form d+I*e, so loop works first time, but second time one of the arguments to the 'Quo' function is of the form d+e _i.

I don't know why you are (apparently?) using two symbols (I and _i) for the same(?) thing, but this looks like a *bad* idea. Without executable code, I cannot investigate this further - and there is a reasonable possibility that I am completely wrong!

@torabi 

was intended to perform very specific simplifications/transformations on a very specific problem in order to reproducea textbook solution (which I failed to do) . Why would you expect these very specific simplifications/transformations to work on a completely different problem?

@Gabriel samaila 

The last two execution groups in the attached will generate exactly the final rows of data in both Table 1 and Table 2 of the reference which you supply. This is good because it suggests that the worksheet is correct.

I have also generated some "multiple" plots altough it is not clear to me which plots you are trying to reproduce - and the supplied reference isn't helpng much!


 

restart

with*plots; ode1 := diff(f(eta), eta, eta, eta)-(diff(theta(eta), eta))*(diff(f(eta), eta, eta))/(theta(eta)-R)-(theta(eta)-R)*f(eta)*(diff(f(eta), eta, eta))/(2*R) = 0

diff(diff(diff(f(eta), eta), eta), eta)-(diff(theta(eta), eta))*(diff(diff(f(eta), eta), eta))/(theta(eta)-R)-(1/2)*(theta(eta)-R)*f(eta)*(diff(diff(f(eta), eta), eta))/R = 0

(1)

ode2 := diff(theta(eta), eta, eta)-N*pr*(diff(f(eta), eta))*theta(eta)+(1/2)*pr*f(eta)*(diff(theta(eta), eta)) = 0

diff(diff(theta(eta), eta), eta)-N*pr*(diff(f(eta), eta))*theta(eta)+(1/2)*pr*f(eta)*(diff(theta(eta), eta)) = 0

(2)

bcs1 := f(0) = 0, (D(f))(0) = 1, (D(f))(11) = 0

f(0) = 0, (D(f))(0) = 1, (D(f))(11) = 0

(3)

fixedparameter := [pr = 100, R = -15.0]

[pr = 100, R = -15.0]

(4)

ode3 := eval(ode1, fixedparameter)

diff(diff(diff(f(eta), eta), eta), eta)-(diff(theta(eta), eta))*(diff(diff(f(eta), eta), eta))/(theta(eta)+15.0)+0.3333333334e-1*(theta(eta)+15.0)*f(eta)*(diff(diff(f(eta), eta), eta)) = 0

(5)

bcs2 := theta(11) = 0, theta(0) = 1

theta(11) = 0, theta(0) = 1

(6)

ode4 := eval(ode2, fixedparameter)

diff(diff(theta(eta), eta), eta)-100*N*(diff(f(eta), eta))*theta(eta)+50*f(eta)*(diff(theta(eta), eta)) = 0

(7)

L := [0, .3, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20]

[0, .3, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20]

(8)

#
# Assorted multiple plots - not really sure which
# plots the OP wants!!
#
  for k to 4 do
      sol_All:= dsolve( eval({bcs1, bcs2, ode3, ode4}, N = L[k]),
                        [f(eta), theta(eta)],
                        numeric,
                        output = listprocedure
                      );
      plt_feta[k]:= plot( (eval(f(eta), sol_All))(eta),
                          eta = 0 .. 1,
                          labels = [eta, f(eta)]
                        );
      plt_fetap[k]:= plot( eval(diff(f(eta), eta), sol_All)(eta),
                           eta = 0 .. 1,
                           labels = [eta, diff(f(eta), eta)]
                         );
      plt_theta[k]:= plot( eval(theta(eta), sol_All)(eta),
                           eta = 0 .. 1,
                           labels = [eta, theta(eta)]
                         );
      plt_thetap[k]:= plot( eval(diff(theta(eta), eta), sol_All)(eta),
                            eta = 0 .. 1,
                            labels = [eta, diff(theta(eta), eta)]
                          ):
  end do:
  plots:-display( [ seq
                    ( plt_feta[j],
                      j = 1 .. 4
                    )
                  ]
                );
  plots:-display( [ seq
                    ( plt_theta[j],
                      j = 1 .. 4
                    )
                  ]
                );
  plots:-display( [ seq
                    ( plt_fetap[j],
                      j = 1 .. 4
                    )
                  ]
                );
  plots:-display( [ seq
                    ( plt_thetap[j],
                      j = 1 .. 4
                    )
                  ]
                );

 

 

 

 

 

#
# Values generated here correspond to those in the
# final row of Table 1 (ie -f''(0) ) in the reference
# supplied by the OP
#
  for k to 4 do
      sol_All := dsolve(eval({bcs1, bcs2, ode3, ode4}, N = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure);
      printf("%16.8f", eval(diff(f(eta),eta,eta), sol_All)(0))
  end do:
 

     -0.47123391     -0.47153126     -0.47192088     -0.47220450

 

#
# Values generated here correspond to those in the
# final row of Table 2 (ie -theta'(0) ) in the
# reference supplied by the OP
#
  for k to 4 do
      sol_All := dsolve(eval({bcs1, bcs2, ode3, ode4}, N = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure);
      printf("%16.8f", -eval(diff(theta(eta),eta), sol_All)(0))
  end do:

      5.54064590      7.59859977     11.16834542     14.92833719

 

 


 

Download someAnswers.mw

@Carl Love 

the OP is employing a simple substitutional code, having more in common with a Caesar cipher than anything vaguely related to the Rivest–Shamir–Adleman algorithm

In the equation you provide, is 'v' actually 'v(t)' ?

The fact that 'v with a dot' appears in the equation implies that 'v' is actually 'v(t)'. But if ir is, then you have a single PDE with two dependent variables y(x,t) and v(t) - generally no way to solve this

be a misprint somewhere in the supplied pdf. The closest I can get to your eqs 3.9.21 and 3.9.22 is shown in the attached.

However notice that (in the attached)

  1. the equation for diff(p(t),t) is of degree 3 in q(t) and degree 2 in p(t)
  2. the equation for diff(q(t),t) is of degree 3 in p(t) and degree 2 in q(t)
  3. the pdf equations are of degree 3 in both p(t) and q(t)

  restart;
  sys:= [ diff(a(t),t)=-mu*a(t)-alpha__6*a(t)*sin(gamma(t))/4,
          a(t)*diff(gamma(t),t)=2*a(t)*sigma-(6/8)*(alpha__1-alpha__2+alpha__3/3)*a(t)^3-alpha__6*a(t)*cos(gamma(t))/2
        ];
#
# Perform the coordinate transformation
#
  trans:=PDEtools:-dchange({a(t)=sqrt(p(t)^2+q(t)^2), gamma(t)=arctan(q(t)/p(t))}, sys, [p(t), q(t)]):
#
# Solve for diff(p(t),t) and diff(q(t),t)
# and do a bit of simplification
#
  ans1:= simplify~
         ( expand~
           ( solve( trans, [ diff(p(t),t), diff(q(t),t)])[] ),
           sqrt,
           symbolic
         );

[diff(a(t), t) = -mu*a(t)-(1/4)*alpha__6*a(t)*sin(gamma(t)), a(t)*(diff(gamma(t), t)) = 2*a(t)*sigma-(3/4)*(alpha__1-alpha__2+(1/3)*alpha__3)*a(t)^3-(1/2)*alpha__6*a(t)*cos(gamma(t))]

 

[diff(p(t), t) = (3/4)*p(t)^2*q(t)*alpha__1-(3/4)*p(t)^2*q(t)*alpha__2+(1/4)*p(t)^2*q(t)*alpha__3+(3/4)*q(t)^3*alpha__1-(3/4)*q(t)^3*alpha__2+(1/4)*q(t)^3*alpha__3-mu*p(t)-2*q(t)*sigma+(1/4)*p(t)*alpha__6*q(t)/(p(t)^2+q(t)^2)^(1/2), diff(q(t), t) = (1/4)*(-3*(p(t)^2+q(t)^2)^(1/2)*(alpha__1-alpha__2+(1/3)*alpha__3)*p(t)^3-2*p(t)^2*alpha__6+8*(p(t)^2+q(t)^2)^(1/2)*((-(3/8)*alpha__1+(3/8)*alpha__2-(1/8)*alpha__3)*q(t)^2+sigma)*p(t)-4*q(t)*mu*(p(t)^2+q(t)^2)^(1/2)-q(t)^2*alpha__6)/(p(t)^2+q(t)^2)^(1/2)]

(1)

#
# Work the above equations individually
# to simplify as much as possible
#
  op(1, ans1[1]) = factor( add( [ op( [2,1..6],ans1[1] ) ] ) )
                   +
                   add( [ op( [2,7..-1], ans1[1] ) ] );
  op(1, ans1[2]) = factor(add([op([2,1..6],expand(ans1[2]))]))
                   +
                   add([op([2,7..8],expand(ans1[2]))])
                   +
                   simplify(add([op([2, 9..10], expand(ans1[2]))]));

diff(p(t), t) = (1/4)*q(t)*(p(t)^2+q(t)^2)*(3*alpha__1-3*alpha__2+alpha__3)-mu*p(t)-2*q(t)*sigma+(1/4)*p(t)*alpha__6*q(t)/(p(t)^2+q(t)^2)^(1/2)

 

diff(q(t), t) = -(1/4)*p(t)*(p(t)^2+q(t)^2)*(3*alpha__1-3*alpha__2+alpha__3)+2*p(t)*sigma-q(t)*mu-(1/4)*alpha__6*(2*p(t)^2+q(t)^2)/(p(t)^2+q(t)^2)^(1/2)

(2)

 

Download simpdiff.mw

@phbmc 

So I suspect that either

  1. It is a version issue: which Maple version are you running?
  2. You are being "mugged" by using 2-D input. Try downloading/running the attached (1-D input) version. If this works, then upload your 2-D input version here using the big green up-arrow in the toolbar, and someone may be able to work out the problem

restart;
L:= proc(c::posint, L::And(list, satisfies(L-> nops(L)>=2*c)))
local p:= L[], A:= p[1+c..-c-1], P:= (p[1..c], p[-c..-1]);
     [seq(
          [P[[seq(p[1..c])]], A, P[[seq(p[-c..-1])]]],
          p= Iterator:-TopologicalSorts(
               2*c, rhs~(op(3, rtable((1..c)$2, (i,j)-> i<c+j, 'storage'= 'triangular'['upper'])))
          )
     )]
end proc:
L(3, [a1,a2,a3, a,a, b1,b2,b3]);

[[a1, a2, a3, a, a, b1, b2, b3], [a1, a2, a3, a, a, b1, b3, b2], [a1, a2, a3, a, a, b3, b1, b2], [a1, a2, a3, a, a, b2, b1, b3], [a1, a2, a3, a, a, b2, b3, b1], [a1, a2, a3, a, a, b3, b2, b1], [a1, a2, b2, a, a, a3, b1, b3], [a1, a2, b2, a, a, a3, b3, b1], [a1, a2, b1, a, a, a3, b2, b3], [a1, a2, b1, a, a, a3, b3, b2], [a1, a2, b1, a, a, b2, a3, b3], [a1, a2, b2, a, a, b1, a3, b3], [a1, b1, a2, a, a, a3, b2, b3], [a1, b1, a2, a, a, a3, b3, b2], [a1, b1, a2, a, a, b2, a3, b3], [a1, a3, a2, a, a, b1, b2, b3], [a1, a3, a2, a, a, b1, b3, b2], [a1, a3, a2, a, a, b3, b1, b2], [a1, a3, a2, a, a, b2, b1, b3], [a1, a3, a2, a, a, b2, b3, b1], [a1, a3, a2, a, a, b3, b2, b1], [a1, a3, b1, a, a, a2, b2, b3], [a1, a3, b1, a, a, a2, b3, b2], [a1, b1, a3, a, a, a2, b2, b3], [a1, b1, a3, a, a, a2, b3, b2], [a3, a1, a2, a, a, b1, b2, b3], [a3, a1, a2, a, a, b1, b3, b2], [a3, a1, a2, a, a, b3, b1, b2], [a3, a1, a2, a, a, b2, b1, b3], [a3, a1, a2, a, a, b2, b3, b1], [a3, a1, a2, a, a, b3, b2, b1], [a3, a1, b1, a, a, a2, b2, b3], [a3, a1, b1, a, a, a2, b3, b2], [a2, a1, a3, a, a, b1, b2, b3], [a2, a1, a3, a, a, b1, b3, b2], [a2, a1, a3, a, a, b3, b1, b2], [a2, a1, a3, a, a, b2, b1, b3], [a2, a1, a3, a, a, b2, b3, b1], [a2, a1, a3, a, a, b3, b2, b1], [a2, a1, b2, a, a, a3, b1, b3], [a2, a1, b2, a, a, a3, b3, b1], [a2, a1, b1, a, a, a3, b2, b3], [a2, a1, b1, a, a, a3, b3, b2], [a2, a1, b1, a, a, b2, a3, b3], [a2, a1, b2, a, a, b1, a3, b3], [a2, a3, a1, a, a, b1, b2, b3], [a2, a3, a1, a, a, b1, b3, b2], [a2, a3, a1, a, a, b3, b1, b2], [a2, a3, a1, a, a, b2, b1, b3], [a2, a3, a1, a, a, b2, b3, b1], [a2, a3, a1, a, a, b3, b2, b1], [a3, a2, a1, a, a, b1, b2, b3], [a3, a2, a1, a, a, b1, b3, b2], [a3, a2, a1, a, a, b3, b1, b2], [a3, a2, a1, a, a, b2, b1, b3], [a3, a2, a1, a, a, b2, b3, b1], [a3, a2, a1, a, a, b3, b2, b1]]

(1)

 

Download carlsCode.mw

@David1

When you set solutions=1, it searches through all the local minima obtained and outputs the best. Hence this setting does not really affect the calculation process, just the output display.  

The above is my interpretation of the following paragraph from the help - my highlighting

solutions = posint -- Set the maximal number of solutions returned. By default the all found solutions are returned. The maximal number of solutions returned cannot exeeds the number of initial points.

Therefore some solutions may "found", but not "returned" Each initial point has the potential to lead to a solution (but might fail):. Furthermore, more than one initial point may lead to the same solution, and such duplicates will not be output

The number of initial points is set according to (my highlighting)

number = posint -- Set the number of initial points simulated. The default value is 100. When initialpoints is provided this option is ignored. The more is number of initial points the more is number of objective function evaluations but also the more are both reliability and accuracy.

So, by default, GlobalSearch will generate 100 initial points and therefore "find" up to 100 solutions: it will "return" all of these found solutions (in order), with the best at the top of the list. If you set solutions=1, it will still find up to 100 solutions, but only "return" 1 - the best.

In the attached worksheet I have tested DirectSearch:-GlobalOptima on a number of  cases from the Wikipedia page

https://en.wikipedia.org/wiki/Test_functions_for_optimization

This page is a useful resource for testing optimizers, although I have no doubt that many other such lists are available, and that different individuals will have their own list of functions for "breaking" optimizers

In each case I have tested both the commands

allSols := DirectSearch:-GlobalOptima(Obj);
oneSol := DirectSearch:-GlobalOptima(Obj, solution=1);

My observations

  1. The 2-D Rastrigin function
    1. allSols[1] agrees with oneSol (within numerical noise) and both agree with the global solution provided by Wikipedia
    2. allSols returns 72 solutions (NB not 100 as expected - presumably some of the initial points led to optima which had aready been found
    3. no significant difference in cpu_usage between allSols() and oneSol()
  2. The Ackley function
    1. allSols[1] agrees with oneSol (within numerical noise) and both agree with the global solution provided by Wikipedia
    2. allSols returns 59 solutions (NB not 100 as expected - presumably some of the initial points led to optima which had aready been found
    3. no significant difference in cpu_usage between allSols() and oneSol()
  3. The 3-D Rosenbrock function
    1. allSols[1] agrees with oneSol (within numerical noise) and both agree with the global solution provided by Wikipedia
    2. allSols returns 1 solution (NB not 100 as expected - suggests that this function is in fact convex, and all initial points leads to the same final optimum)
    3. no significant difference in cpu_usage between allSols() and oneSol()
  4. The Beale function
    1. allSols[1] agrees with oneSol (within numerical noise) and both agree with the global solution provided by Wikipedia, although in both cases, there is a warning that more function evaluations may be necessary
    2. allSols returns 3 solutions (NB not 100 as expected - presumably there are only three minima and all initial points lead to one of these three)
    3. no significant difference in cpu_usage between allSols() and oneSol()
  5. The GoldStein-Price function
    1. allSols[1] agrees with oneSol (within numerical noise) and both agree with the global solution provided by Wikipedia
    2. allSols returns 3 solutions (NB not 100 as expected - presumably there are only three minima and all intial points lead to one of these three)
    3. no significant difference in cpu_usage between allSols() and oneSol()

Worksheet is saved with output which makes a bit lengthy to display inline here, hence just the link

testOptim.mw

@David1 

of the 'solutions=1' option - and I admit that it is dfficult to justify from help pages, but bear with me

The way the GlobalSearch() command works is to start form a *lot* of different initialValues: from each initialValue, keep iterating until a minimum is found; technically this will be a local minimum. Having started from a *lot* of different initiialValues, a *lot* of (possibly) different local minima will be obtained. When you set solutions=1, it searches through all the local minima obtained and outputs the best. Hence this setting does not really affect the calculation process, just the output display.

Is it a definitely a global minima? No!

Does any setting of the 'solutions' option guarantee a global minimum? No!

Because as I keep repeating there is no known way to guarantee a global minima for anything other than a convex function

post code using the big green up-arrow in the MaplePrimes toolbar. No-one is going to retype your code from a 'picture' - it's error-prone and too time-consuming.

Idle curiosity, what happens if you add the option 'numeric' to the Quantile() command?

@waseem 

Apart from anything else, you have introduced "square brackets" into your desired form - in Maple these represent indexable quantities. So your desired form now contains several, one-element, lists - whihc you probably don't want!

In general, geting exactly the output which you think is 'prettiest' or 'most desirable' from Maple can take more time than it is worth. Having said that, the attached gets 'close' to the form you want

restart;
u[1](r,z):=(C[o]^2*exp(2*lambda*z)-D[o]^2*exp(-2*lambda*z))+(lambda/2)*(C[o]*exp(lambda*z)-D[o]*exp(-lambda*z))*(1-r^4)+A[1](z)*r^2;
M1:=expand(int(u[1](r,z), z));
MM1:=add(simplify~([op([1..5],collect(simplify(M1, power), [r, lambda]) )]));

C[o]^2*exp(2*lambda*z)-D[o]^2*exp(-2*lambda*z)+(1/2)*lambda*(C[o]*exp(lambda*z)-D[o]*exp(-lambda*z))*(-r^4+1)+A[1](z)*r^2

 

(1/2)*C[o]^2*(exp(lambda*z))^2/lambda+(1/2)*D[o]^2/(lambda*(exp(lambda*z))^2)-(1/2)*r^4*C[o]*exp(lambda*z)-(1/2)*r^4*D[o]/exp(lambda*z)+(1/2)*C[o]*exp(lambda*z)+(1/2)*D[o]/exp(lambda*z)+r^2*(int(A[1](z), z))

 

-(1/2)*r^4*(C[o]*exp(lambda*z)+D[o]*exp(-lambda*z))+r^2*(int(A[1](z), z))+(1/2)*C[o]*exp(lambda*z)+(1/2)*D[o]*exp(-lambda*z)+(1/2)*(C[o]^2*exp(2*lambda*z)+D[o]^2*exp(-2*lambda*z))/lambda

(1)

 


 

Download intProb2.mw

@mmcdara 

Starting with 1978 as base year and counting by two's the five year overage global temperature. where 

 

the temperature is given by the data : x[0],...,x[12] 

so x[0] corresponds to 1978 which makes x[30] correspond to 2038

@mmcdara 

Most of your contribution in the reply

TimeSeries or not TimeSeries?

and your reference to your earlier  TimeSeries.mw is invalidated by the simple fact that the OP's data is in 2-year steps so the quantity x[30] refers to 60 years after the start data of 1978, ie 2038. So why does the analysis in TimeSeries.mw stop at 2008 only 30 years after the start date?? OP's original question states quite clearly " Starting with 1978 as base year and counting by two's"

If you check my last post (where I use the time series approach, you will note that all data is presented with two year intervals (as have all my earlier posts)

If you want to correct your all of your previous analyses to account  for this simple fact, your answers will get very close to those I have already given

The simple reason I don't propose to add much further to this thread is that I am not a great believer in lengthy extrapolation from actual data. Such extrapolation is always mathematically possible - but tiny discrepancies in fitted functions can lead to huge discrepancies in extrapolated values, so is the process meaningful/worth it

@mmcdara 

about necessarily doing "more reliable" things with the Timeseries() package. I too have never used it, but I gave it a quick try in the attached - whilst the results "look plausible", the required extrapolated values come out noticeably higher than those I obtained previously.

I haven't got the time/energy/expertise/patience to examine more fully. If I were doing a lot of statistical analyes (which I'm not!) then it might be a worthwhile exercise to investigate more fully

oddFit3.mw

First 81 82 83 84 85 86 87 Last Page 83 of 207