thx for nailing it down, that was worth the effort

May be that comes through a bug in 'simplify' or 'abs'.

Define I3:=PDEtools[dchange](x=1/u ,I1) with integrand p:=op(1,%) =
abs(-1/u*sin(u^2)+cos(u^2)*u) / u^2, which is your I2, but with the
absolute sign is not forgotten.

Now apply 'simplify' to I3:

  simplify(I3);

                    infinity
                   /               2         2   2
                  |          -sin(u ) + cos(u ) u
                  |          --------------------- du
                  |                    3
                 /                    u
                   1

which Maple evaluates to -1/2*sin(1) and that coincides with the
given value for I0 ( u=1/x shows it to - I0).



This leads to find the following bug: 

Both simplify(p) and simplify(p,symbolic) behave correct. However
they do not if using assumptions:

  q;


                      |        2              |
                      |   sin(u )        2    |
                      | - ------- + cos(u ) u |
                      |      u                |
                      -------------------------
                                  2
                                 u

  simplify(q) assuming 0<u;
  # plot(%,u=1 .. sqrt(2*Pi));

                              2         2   2
                        -sin(u ) + cos(u ) u
                        ---------------------
                                  3
                                 u

To be sure that this is not caused through the 'assuming' command
one produces the same error through 'assume':

  assume(0<v): eval(p, u=v): simplify(%);

                              2         2   2
                        -sin(v ) + cos(v ) v
                        ---------------------
                                  3
                                 v

Strange that eval(p, u=sqrt(v)), eval(p, u=v^2) or eval(p, u=(v+1))
work as expected.

Since it is not unlikely that IntegrationTools[Change] uses 'simplify'
that would 'explain' it and would not mean that the integration itself
has an error here: may be 'Change' wants to use information about the
range of integration.


In fact the reason might be in the 'abs' function:

  r:=(sin(v^2)+cos(v^2)*v^2)/v^3; # we have 0 < v
  s:=abs(r);
                                2         2   2
                           sin(v ) + cos(v ) v
                      r := --------------------
                                     3
                                    v


                                2         2   2
                           sin(v ) + cos(v ) v
                      s := --------------------
                                     3
                                    v

Which is false, just plot it. While eval(r, v=u); abs(%); works ok ...

only a minor comment:

  PDEtools[dchange](x=1/u ,I1);
  student[changevar](x=1/u ,I1,u);

seem not to ignore the absolute value.

Have not used VBA for a longer time ... First correct a 'feature' in Excel's
modulo function and then translate the EGCD given by Roman Pearce. The code
has just to be copied into a module of an Excel VBA project.


Function myMod(a As Long, p As Long) As Long
' check that for a=-20, p=7: Excel would return -6, we want positive numbers
myMod = (a) Mod p
If a < 0 Then
  myMod = myMod + Abs(p)
End If
End Function


Function invMod(m As Long, p As Long)
' returns invers modulo p using the extended GCD algorithm (Knuth?)
Dim a As Long, b As Long, q As Long, t As Long, x As Long, y As Long

a = p
b = m
x = 1
y = 0
While Not (b = 0)
  t = b
  q = Application.WorksheetFunction.Floor(CDbl(a / t), 1)
  ' CodeGeneration[VisualBasic] is not aware of that VBA, it uses VB
  b = a - q * t
  a = t
  t = x
  x = y - q * t
  y = t
Wend
If (y < 0) Then
  y = y + p
End If
invMod = y
End Function


Sub tst_invMod()
' test, to be displayed in the debug window
Dim a As Long, p As Long, b As Long, ab As Long, j As Long

p = 23 ' prime
For j = 1 To p - 1
  a = j
  b = invMod(a, p)
  ab = myMod(a * b, p)
  Debug.Print a & "^(-1) mod(" & p; ") = " & b, _
    "a*b mod(" & p & ") = " & myMod(a * b, p)
Next j
End Sub


The largest (signed) integer of type long is 2^31 in VBA. Taking
p = prevprime(2^31) = 2147483647, a = prevprime(p) = 2147483629
one immediately gets a^(-1) = 119304647, which equals Maple's answer.

The ugly worksheet & search method is limited to p < 65536 and
should be quite slow. 

but a fixed correlation does not make sense, it is known to by stochastic & clustered

portfolio simualtion is a field of its own, especially if one has assets of different types (say credit + equity)

it is minor a question of Maple or not ...

here some reading www.nuclearphynance.com/User%20Files/2/Dispersion%20-%20A%20guide%20for%20the%20clueless%201.1.pdf

and as a side remark: the subprime crisis shows how much it is worth to believe in models and simulation :-)

it means that one has to sum over the roots of the equation, which is a polynomial of degree 3 - so they are 'explicitly' known

and usually (depending on your setting) will supress those complicated answers, you can enforce it to show them by 'allvalues'

so may use allvalues(%):  evala(%):  simplify(%,size); to see the result

I played with it, using Digits:=14, here is my hack:

  F:=lhs(eq1)-rhs(eq1); 
  RootOf(F,lambda, 1e-10 .. 10 ); 
  R:=unapply(%,alpha);

Then plot(R(alpha), alpha=0..6) has a some 'jumping' just below 5
and in other plot ranges.

Looking at the Taylor series by MultiSeries:-series(%, lambda, 3)
lets me try to kick off the solution lambda=0 by re-defining

  G:=F/lambda^(1/2*alpha-1/2) # same zeros as F

Proceeding similar for G still gives a similar problem. So try better:

Now up to 1st Order MultiSeries:-series(G, lambda, 2) = 0 is solved
by alpha. And one can guess that there are 2 regimes: lambda ~ alpha
below 2 and lambda ~ alpha/3 for large values.

Thus I try:

  S:=proc(a)
    global G, lambda;
    eval(G, alpha=a);
    fsolve(%, lambda=a/4);
  end proc;

Then some test

  S(2); 
                           2.5064979940915
  S(25.3);
                           9.5380145635235

and the following plots smoothly:

  plot('S'(a), a=0..20);

what does it contain? i do not want to download huge things + install + clicking through topics ...

well, I would not want to test it on my own PC, but if i have some time may be i will try to install and test it at work (this usually means to buy some icecream for the admins and otherwise keep them busy by complaining formatting problems while printing letters on various printers through the company or so ...)

in former times i used to change system times for testing several Excel solutions (to see, whether they work correctly beyond certain date switches having a complicated trading calendar) or old sheets (having time values to be set only through system time, back from a productive environment) or Word macros (dito) or compiled DLLs (similar).

That would totaly mess up any Maple installation , no?

I prefer to have Sum and use lhs(%) - rhs(%) with is(%=0):  

  y*(Sum(a[n]*y^(n-2)*n*(n-1), n = 0 .. infinity)) 
  - 
  (Sum(a[n]*y^(n-1)*n*(n-1), n = 0 .. infinity));
  combine(%);
  simplify(%);
                                  0

For the database: twice continuously differentiable may be needed to work *properly* (but have not looked it up).
At least you can try to approximate your function by such things?

The other is more work and related to the (by numerical errors only?) above example: you could try to sum the
squares of the components of your equations. And try to use a minimizer, since a minimum =0 gives a zero,
the minimzer should allow to search not only locally. But I am to lazy to do that.

I was just playing with it again, for P = 0 there are various solutions.

  eval(eqs,q=0); fsolve(%, indets(%,symbol)); 
  fnormal(%): simplify(%,zero); # this cleans up the numerical results
  R:=identify(%);               # and gives a nice appearance

     {C19 = 0, C4 = C4, C5 = -C8, C8 = C8, C13 = -C16, C12 = C12,
      C16 = C16, C20 = C20, C24 = C24, C28 = C28, C32 = C32, C7 = 0,
      C31 = 0, C29 = -C32, C15 = 0, C23 = 0, C17 = -C20, C25 = -C28,
      C11 = 0, C3 = 0, C27 = 0, C1 = -C4, C9 = -C12, C21 = -C24}

  eqs: eval(%,q=0); eval(%, R); # shows the system is satisfied

  equations: eval(%,P=0); eval(%, R); # dito, for original problem

The solution show some stuff like "C4=C4". This usually means, that those
variables can be choosen independently.

Let us try that through 2 examples:

Set C4=2 in R, feed it into your system and use P=0. Then all except the
C4 is replaced. Now solve for that (and see it gives C4=2 as desired):

  eval(R, C4=2): eval(equations, %): eval(%, P=0); 

     {-2. + C4 = 0, -0. = -0., 0 = 0., 0 = -0., 0. = 10. - 5 C4,
      -2. + C4 = 0., -0. = 0, 0. = 0, 0. = 0., 0 = 0}

  solve(%, C4);
                              {C4 = 2.}

  # just for fun ...
  eval(R, [C4=Pi, C8=sqrt(2)]): eval(equations, %): eval(%, P=0);  
                      1/2            1/2
     {-1. Pi + C4 = -1. 2    + C8, -1. 2    + C8 = 0., 0. = 5. Pi - 5 C4,

      -0. = -0., 0 = 0., -1. Pi + C4 = 0, 0 = -0.,
                1/2
       0. = 5. 2    - 5 C8, 0. = 0, 0. = 0., -0. = 0, 0 = 0}


  solve(%, {C4,C8}); identify(%);

      {C4 = 3.1415926535898, C8 = 1.4142135623731}

                                1/2
                         {C8 = 2   , C4 = Pi}

A check shows that this is (of course?) correct, using

  eval(R, [C4=Pi, C8=sqrt(2)]): eval(equations, %): eval(%,tmp): eval(%, P=0);

     {-0. = -0., 0 = 0., 0 = -0., 0. = 0, 0. = 0., -0. = 0, 0 = 0}



Now ... I think that especially menas, that you will have problems to judge
the situation for P = sqrt(q) close to zero, at least numerical it is odd:
Choose some q (=dummy) of small size, insert it into eqs, call it tmp:

  dummy:=1/Float(1, floor(Digits*0.4) + 1);
  'eval(eqs,q=dummy)'; 
  tmp:=%;  
  #indets(tmp,symbol); nops(%); nops(tmp);

Now we lost 1 indeterminate, but fsolve wants #equations = #indeterminates.
Staring at the system lets me try to remove the first one and to solve for
the rest - which are the C's, the solution set is named M:

  tmp: convert(%,list): {seq( %[i], i=2 .. nops(%))}: 
  M:=fsolve(%, indets(%,symbol),fulldigits); lprint(%);

Then I feed it into the equations, set q = dummy to see it is satisfied
(writing the equations as lhs - rhs = 0), at least numerical (as allays here):

  'eval(equations, M)'; 'eval(%,P=dummy^2)'; evalf(%): 
  map( t -> lhs(t)-rhs(t),%); 
  fnormal(%);
  
                    eval(equations, M)|
                                      |         2
                                      |P = dummy

                              {0., -0.}

Doing the same with P=dummy^2/2 this is *not* seen as solution. Which says
that it is not simply satisfied, because P is small.

However using 0 .. 1e-4 as range for q in fsolve for a direct result did not
terminate within a time I was willing to wait.

So I think you need an additional reason, why the above practical does not
apply, if you want to use P ~ 1.66 as appropriate value.

Using the differential operator 'D' (may be using convert(expr,D) first)
it is simply D(u*v).

Explanation: Maple seems to 'understand' how basic operations have to be
applied to functions.

  t:='t':

  someConstants:=[0,1,2,evalf(Pi), Pi];
  for b in someConstants do
    p:=b;
    print('p'=p, 'p(t)'=p(t));
  end do:

                           p = 0, p(t) = 0

                           p = 1, p(t) = 1

                           p = 2, p(t) = 2

                 p = 3.141592654, p(t) = 3.141592654

                         p = Pi, p(t) = Pi(t)

So if the constant is not a symbol like Pi an assignment like p:=2 is 
understood as a function taking that constant as value. And no, it would
not work for the Natural for example:

  assume(n::integer): getassumptions(n);
  p:=n;
  p(t);

                                p := n

                                 n(t)

The identity one has to define extra by (id: x -> x works very reliable):

  id:=unapply(x,x);

                             id := x -> x

  'id( t )': '%'=%;
  'id( exp(sin(tau)) )': '%'=%;

                              id(t) = t

                  id(exp(sin(tau))) = exp(sin(tau))


This also works (formally) for the basic arithmetic operators:

  ArithmeticOperators:=[`+`, `-`, `*`, `/`, `^`];
  for b in ArithmeticOperators do
    p:=b(u,v);
    print('p'=p, 'p(t)'=p(t));
  end do:

                    p = u + v, p(t) = u(t) + v(t)

                    p = u - v, p(t) = u(t) - v(t)

                      p = u v, p(t) = u(t) v(t)

                                         u(t)
                         p = u/v, p(t) = ----
                                         v(t)

                            v             v(t)
                       p = u , p(t) = u(t)


Thus u*v is understood as the function t -> u(t)*v(t) and since Maple
expand automatically (otherwise: enforce it) one can simply write:

  D(u*v);
  %(t);
  convert(%,Diff);
                           D(u) v + u D(v)

                     D(u)(t) v(t) + u(t) D(v)(t)

                   /d      \             /d      \
                   |-- u(t)| v(t) + u(t) |-- v(t)|
                   \dt     /             \dt     /


The last command is onlyto show it in more common notation.

I think it does not, since it picks just one ...

Here is what I did (within your sheet using Digits:=14 [which has
no extra cost usually])):

  # the equations are in sqrt(P), eliminate that
  equations:
  convert(%,rational):
  eval(%,P=q^2): simplify(%) assuming 0 < q: evalf(%):
  
  eqs:=%; # now sqrt(P) is just q, we want to know q


  eqs:
  fsolve( %, indets(%,symbol), q=0..sqrt(10));
  select(has,%,q);
  rhs(op(%))^2; # P = q^2

                        {q = 1.2893925628043}

                           1.6625331810150

This seems to be returned quite quickly (which is your P).


For iterations (may be not needed here) I would try something like

  eqs:
  fsolve( %, indets(%,symbol), q=0..sqrt(10));
  select(has,%,q);
  rhs(op(%)):
  P=%^2;
  qNew:=%%;
  
  # use it as new starting value and a limit for
  # computational time
  sols:=timelimit(10,
    fsolve( eqs, indets(eqs,symbol), q=0..qNew));
  
  select(has,sols,q);
  rhs(op(%)):
  P=%^2;
  qNew:=%%;


Can not see, whether your system has a (constant?) periodics in P.


I think the general problem is (even for your case over the Reals):

If #equations = #variables, then the generic solution constists of
isolated points. One can neither expect that these are finite in
number (say: you have periodics). Nor does one know a priori that
all are catched even for the compact (=finite) case.

And an additional condition may easily mean: you have to search the
space of solutions (=closest to zero and positive, real).

If your system is algebraic (=given by polynomial equations) then
you might try the Gröbner basis methods in Maple, to describe the
zero dimansional variety (no, I never seriously worked with it).

Since your system is analytical for that you either would transform
it to an algebraic one (if ever possible) or approximate it (not so
claer what has to be done, except you know some location for your
desired solution), do not know of Maple has s.th. for analytical
given zero-dimensional varieties (=discrete points as intersection
of analytical functions) - even if that is mathematical equivalent
to an algebraic situation (locally).


For short without pseudo-smart blabla: seems you have to search the
solutions for the appropriate one, no guarantee for completeness.