Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 181 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

Corrected code...

tomleslie 13876
July 11 2018
0 0

Fixing the problem described by Kitonum, and generally tidying up the code you presented results in the code/graphs shown in the attached.

You claim you want to "vary more than one variable at a time". Which variables do you want to vary?

NULL

> 

restart; w := .572433; x := 1.32156; y := 5.29387; z := .853043; b := -.1; pr := 6.2; m := .5; N := [10, 20, 30]; K := [black, red, green, blue, green, orange, gold, gray]; L := [solid, solid, solid, dash, solid, dot, dash]; odes := diff(f(eta), `$`(eta, 3))+w*x*(f(eta)*(diff(f(eta), `$`(eta, 2)))-(m*m)*(diff(f(eta), eta))-(diff(f(eta), eta))^2) = 0, y*(diff(theta(eta), `$`(eta, 2)))/(pr*z)-b*f(eta)*(diff(f(eta), eta))*(diff(theta(eta), eta))-b*f(eta)^2*(diff(theta(eta), `$`(eta, 2)))+f(eta)*(diff(theta(eta), eta)) = 0; bcs := (D(f))(0) = 1, (D(f))(20) = 0, theta(0) = 1, theta(20) = 0

diff(diff(diff(f(eta), eta), eta), eta)+.7565045555*f(eta)*(diff(diff(f(eta), eta), eta))-.1891261389*(diff(f(eta), eta))-.7565045555*(diff(f(eta), eta))^2 = 0, 1.000946025*(diff(diff(theta(eta), eta), eta))+.1*f(eta)*(diff(f(eta), eta))*(diff(theta(eta), eta))+.1*f(eta)^2*(diff(diff(theta(eta), eta), eta))+f(eta)*(diff(theta(eta), eta)) = 0

 (1)
> 

for j to nops(N) do sol1 := dsolve([odes, bcs, f(0) = N[j]], numeric, maxmesh = 1024); fplt[j] := plots[odeplot](sol1, [eta, f(eta)], color = K[j], axes = boxed, linestyle = L[j]); tplt[j] := plots[odeplot](sol1, [[eta, theta(eta)]], color = K[j], axes = box, linestyle = L[j]) end do; plots:-display([seq(fplt[j], j = 1 .. nops(N))]); plots:-display([seq(tplt[j], j = 1 .. nops(N))])

 

 
> 

 

NULL

Download odeProb.mw

 

Not sure what to do next...

tomleslie 13876
July 11 2018
0 0

@Bosco Emmanuel 

If I run this latest code on the oldest version of Maple which I have (Maple 18, released 2014, approximately), then I too get the message about "no iterations performed because....". This does not happen on the latest (ie current) version, Maple 2018, where as I stated before, the Minimize() command runs for a while before I get bored and then stop it.

Does Maple 18 *think* that the initial point is (at least a local) minimum - maybe? A simple test is just to change Minimize() to Maximize() and see what happens. From the same initial point, this produces a "maximum" of 305.535138631842926  at the point 2.99999999999855, 1.00000063538814, 0.500000000007119, 0.299999159777189, 1.00000037825258, 1.50000000000641

This suggests trying to Maximise() the reciprocal of the original function - don't bother: it returns the same error message "no iterations performed because....".

So I did a little more investigation, with a simple printf() statement. With the original code, Minimize() calls the function 'lc' 117 times, for its "first" iteration, with slightly different values of the parameters, more or less as I would expect. On each of these function calls, three new set of random variables will be generated, even when the relevant parameter p1, p2, p3 doesn't change. So each function call will return a new value for 'Error', depending on both the particular set of random variables which is generated, and the the values of the parameters. You can prove this to yourself just by checking multiple evaluations of 'lc()' with the same set of parameters, as in

lc(3, 1, .5, .3, 1, 1.5);
lc(3, 1, .5, .3, 1, 1.5);
lc(3, 1, .5, .3, 1, 1.5);
lc(3, 1, .5, .3, 1, 1.5);

which on my machine returns  253.8115879115319, 255.98597940362566, 255.48304268485603) 281.5371300272584. I can see why this would lead to all sorts of "issues" for an optimisation routine, because the effect of changing one parameter in a certain direction might be reversed just by the 'new' random variable set.

Just to double-check how serious this behavior might be, I ran the same DirectSearch() command four times with an identical initial point. Results were

 Err:= 37.4781849722149  pars := [5.03878871903693, 0.0143456613127857,  0.0571605855762387, 0.407489231463475, 0.974677058139177 0.960498467662941]
 Err := 2.84343381535019 pars:= [4.66659603431891, 43.1941143025958, 1.38203637080510, 0.584475219434157, 0.627722132098533, 0.619158462783508]
 Err := 2.48150382650395 pars := [4.66920238506976, 1445.07011203137, 1.12331530564846,  0.717627279664310, 0.644897160956551, 1.39303150886333]
 Err := 36.3338681532953 pars := [4.98661618524935, 0.00531500177051038, 0.0237137981290860, 0.417693660380795, 2.29048947849934, 0.401341881397136]

Notice the huge variation in the "minimum" error obtained, and also the variation in the parameters which lead to roughly similar minima.

The only conclusion I can draw is that this "optimisation" is a complete waste of time (and changing the optimiser won't help).

I return to my previous observation about an explanation of the precise nature of your problem, and exactly what you are trying to achieve - in maths/English, not code, because whatever the problem actually is, you chances of solving it with your preent approach is precisely nil.

 

What are you trying to achieve????...

tomleslie 13876
July 10 2018
0 0

@Bosco Emmanuel 

Please read this very carefully - bedore you do anything with the attached code

  1. I didn't understand the logic behind your calculation before, and I understand it even less as a "function-to-be-minimized". Forget about Maple, coding etc. Can you describe what you are attempting in mathematical terms? Because whatever you are trying to achieve, I am more and more certain that you are going about it the "wrong" way.
  2. Can I make your latest worksheet run? Yes (you introduced a few syntax errors with your attempt to use the Optimization package, but these I can fix, see later). However the Minimize() command doesn't terminate in the time I was prepared to wait (about 10mins). I wasn't really surprised. Maple's built-in optimisation routines are not particularly robust. Which is probably why Maple can get away with selling an add-on package called the Global Optimization Toolbox. I'm assuming you don't have this add-on, I don't because I can't jusify the expense.
  3. There is another free add-on Optimization package for Maple called DirectSearch, which is generally considered to be a lot more robust than the built-in Optimization routines.  This is available from the Maple Application centre, so has to be downloaded/installed. I have no idea whether this package will run in Maple12 (it was developed using Maple13, and has been subsequently modified, so I doubt it). So before you even think about this option - please consider my remarks in (1) above
  4. I decided to try your problem using the DirectSearch add-on. You will not be able to run this because you won't have the DirectSearch() package installed. This produced an answer in about 10secs- see the attached
    1. With no "Optimization", Running the routine lc() with your initial point [3, 1.0, .5, .3, 1, 1.5] returns Error := 253.811587909557
    2. The DirectSearch() package returns Error := 37.4781849722149, at the point [5.03878871903693, 0.143456613127857e-1, 0.571605855762387e-1, .407489231463475, .974677058139177, .960498467662941]
    3. Since your specified optimisation process requires the generation of random numbers, I am doubtful about how reproducible this result will be - again I think that your problem is ill-posed
  5. If you want to try the built-in Minimize() for longer, then just uncomment it in the attached (and comment out the subsequent DirectSearch commands), Just be ready to terminate it manually!
  6. There are several "inefficiencies" in your code which I have not attempted to correct. Since optimisation is an iterative process, the procedure lc() has to be executed multiple times. As written, this means that the input data file will be 'read' multiple times. This is time-consuming, and completely pointless, but I haven't bothered to fix this (or other inefficiencies) for fear of introducing something else you will be unable to run
> 

restart; with(Statistics); with(Optimization); lc := proc (p1, p2, p3, p4, p5, p6) local RT, TC, N, T, z, n, l, B, p; p := [p1, p2, p3, p4, p5, p6]; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1.0/p[1]), 1001), Sample(Exponential(1.0/p[2]), 1001), Sample(Exponential(1.0/p[3]), 1001)], scan = columns); z := ImportMatrix("C:/Users/TomLeslie/Desktop/Tom.txt", source = delimited, delimiter = "\t", format = rectangular, datatype = sfloat, transpose = false, skiplines = 0); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Error := lc(3, 1, .5, .3, 1, 1.5); ans := DirectSearch:-Search(lc, initialpoint = [3, 1.0, .5, .3, 1, 1.5], assume = positive); Error := ans[1]; parVals := convert(ans[2], list); lc(parVals[])

HFloat(253.81158790955718)

  

HFloat(37.47818497221495)

  

[HFloat(5.038788719036934), HFloat(0.014345661312785718), HFloat(0.05716058557623869), HFloat(0.4074892314634753), HFloat(0.9746770581391767), HFloat(0.9604984676629409)]

  

HFloat(48.59838971582652)

 (1)
> 

NULL

Download ds1.mw

You mean...

tomleslie 13876
July 09 2018
0 0

@waseem 

something more like the plot shown in the attached

plotProb2.mw

I just...

tomleslie 13876
July 09 2018
0 0

@Bosco Emmanuel 

downloaded it from the link in my previous post (exactly as you would do) - and it downloaded and executed as normal. Can't see why you would have a problem.

However I will upload again


 

> 

  restart;
  with(Statistics):
#
# Initialise a few things
#
  p:= [3, 1, 0.5, 0.3, 1, 1.5]:
  RT:= Array(1..3, p[4..6]):
  TC:= Array(0..1001, fill=0):
  N:= Array(1..3, [0, 0, 0]):
  T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001),
                         Sample(Exponential(1/p[2]), 1001),
                         Sample(Exponential(1/p[3]), 1001)
                       ],
              scan=columns
            ):
#
# Do the calculation
#
  for n to 1001 do
     # l:= min[index](T[..,n]):
      if   T[1,n]<=T[2,n] and T[1,n]<=T[3,n]
      then l:=1
      elif T[2,n]<=T[1,n] and T[2,n]<=T[3,n]
      then l:=2
      else l:=3
      fi:
      N[l]:= N[l]+1;
      B[N[l], l]:= T[l,n]+RT[l]; # changed T[1,n] to T[l,n]
      TC[n]:= TC[n-1]+B[N[l], l]
   end do:
   z:= ImportMatrix("C:/Users/TomLeslie/DeskTop/Tom.txt"):
#
# Get the result
#
   #Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000);
    Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000);

HFloat(253.8115879115319)

 (1)
> 

###############################################################
#
# Define procedure which does the same
# calculation as the above
#
  restart:
  lc:=proc( p::list )
          #
          # Declare/initialise local variables
          #
            local RT:= Array(1..3, p[4..6]),
                  TC:= Array(0..1001, fill=0),
                  N:= Array(1..3, [0, 0, 0]),
                  T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001),
                                         Sample(Exponential(1/p[2]), 1001),
                                         Sample(Exponential(1/p[3]), 1001)
                                       ],
                              scan=columns
                            ),
                  z:= ImportMatrix("C:/Users/TomLeslie/DeskTop/Tom.txt"),
                  n, l, B;
          #
          # Ensure Statistics package is loaded
          #
            uses Statistics:
          #
          # Do calculations
          #
            for n to 1001 do
              #  l:= min[index](T[..,n]):
                if   T[1,n]<=T[2,n] and T[1,n]<=T[3,n]
                then l:=1
                elif T[2,n]<=T[1,n] and T[2,n]<=T[3,n]
                then l:=2
                else l:=3
                fi:
                N[l]:= N[l]+1;
                B[N[l], l]:= T[l, n]+RT[l]; # changed T[1,n] to T[l,n]
                TC[n]:= TC[n-1]+B[N[l], l]
            end do:
          #
          # Return result
          #
            return add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000);
       end proc:
#
# Run the calculation with the same data
# as used in the above
#
  Error:=lc([3, 1, 0.5, 0.3, 1, 1.5]);

HFloat(253.8115879115319)

 (2)
> 

############################################
# A 'Quirk' of sample generation
#
# Note that in the following two methods
# the 'first sample always agrees, but
# subsequent samples do not
#
  restart;
  with(Statistics):
  Sample(Exponential(1/0.3), 5);
#
  restart;
  with(Statistics):
  X:=Exponential(1/0.3):
  for i from 1 by 1 to 5 do
      Sample(X,1)[1];
  od;

Vector[row](%id = 18446744074335814230)

HFloat(13.726161926962995)

HFloat(4.783158197400881)

HFloat(1.510645808521173)

HFloat(11.110290254440143)

HFloat(1.5214776243071664)

 (3)
> 

restart; with(Statistics); p := [3, 1, .5, .3, 1, 1.5]; RT := Array(1 .. 3, p[4 .. 6]); TC[0] := 0.; N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001)], scan = columns); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; z := ImportMatrix("C:/Users/TomLeslie/Desktop/Tom.txt"); Error := add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000)

HFloat(253.8115879115319)

 (4)
> 

restart; with(Statistics); lc := proc (p::list) local RT, TC, N, T, z, n, l, B; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Statistics:-Sample(Exponential(1/p[1]), 1001), Statistics:-Sample(Exponential(1/p[2]), 1001), Statistics:-Sample(Exponential(1/p[3]), 1001)], scan = columns); z := ImportMatrix("C:/Users/TomLeslie/Desktop/Tom.txt"); for n to 1001 do if T[1, n] <= T[2, n] and T[1, n] <= T[3, n] then l := 1 elif T[2, n] <= T[1, n] and T[2, n] <= T[3, n] then l := 2 else l := 3 end if; N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Error := lc([3, 1, .5, .3, 1, 1.5])

HFloat(253.8115879115319)

 (5)
> 

 


 

Download randStuff4.mws

More extensive testing...

tomleslie 13876
July 08 2018
0 0

@Carl Love 

I have just tried both of the worksheets Kitonums.mw and mine.mw in all of the Maple versions I have with the following results

Maple Version     Kitnonums     MIne

         18                  works        works
      2015.2              works        works
      2016.2              works        works
      2017.3              works        works
      2018.1              Errors       works

In my original post, I only tried Maple 2018.1 - hence my observation about  Kitonums.mw producing errors, and mine.mw not producing errors. Maybe I should have checked these with earlier Maple versions, but I didn't: mea culpa, etc

For my earlier post, whilst trying to debug the failure of Kitonums.mw in Maple 2018.1, I convinced myself that it was something to do with the use of midline 'dot' operators for matrix multiplication. I cannot now remember why I came to this conclusion, since this syntax *ought* to be OK.

However if you check the differences between these two worksheets, AFAIK the only significant discrepancy is the way matrix/vector/matrix multiplications are handled. The worksheet mine.mw uses the incredibly long-winded MatrixMatrixMultiply(), VectorMatrixMultiply() and MatrixVectorMultiply() calls, rather than a midline 'dot'. Is this why it works? I have no idea!

Just for clarity, all of the above was performed on 64-bit Windows 7



 

Not for me it doesn't...

tomleslie 13876
July 08 2018
0 0

@Kitonum 

If I execute this worksheet in Maple 2018.1, it throws a couple of errors :-(

See the attached


 

> 

restart; with(LinearAlgebra); kernelopts(version); printlevel := 2; M := 3; k := 1; MFD := 2; `&vartheta;0` := 4/3; `&vartheta;1` := 3/4; `&vartheta;2` := 3/2; Delta := 2^(k-1)*M; nu := 1; for i while i <= 2^(k-1) do for j from 0 while j <= M-1 do `&psi;&psi;`[M*i+j-M+1] := simplify(2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+nu)*GAMMA(nu)^2/(Pi*2^(1-2*nu)*GAMMA(j+2*nu)))*(sum((-1)^i1*GAMMA(j-i1+nu)*(2*(2^k*x-2*i+1))^(j-2*i1)/(GAMMA(nu)*factorial(i1)*factorial(j-2*i1)), i1 = 0 .. floor((1/2)*j)))) end do end do

`Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769`

 (1)
> 

NULL
> 

r := 2^(k-1)*M; print("Gegenbauer wavelets vector having order" = 2^(k-1)*M); `&psi;x` := `assuming`([simplify(combine(Vector[column](r, proc (i) options operator, arrow; `&psi;&psi;`[i] end proc)), radical)], [nu > 0]); `&psi;t` := eval(`&psi;x`, x = t); print("Unkown vector need to determine with order" = 2^(k-1)*M); C := Matrix(r, r, symbol = a); print("Trial solution of GWM when k" = k, "and M" = M); u := `assuming`([simplify(combine(Transpose(`&psi;x`).C.`&psi;t`), radical)], [nu > 0]); for j1 while j1 <= MFD do B[j1] := proc (M) Matrix(M, M, {seq((i, i) = factorial(i-1)/GAMMA(i-`&vartheta;` || j1), i = j1+1 .. M)}) end proc end do; for j1 while j1 <= MFD do FD[j1] := B[j1](M) end do; for j1 while j1 <= MFD do B1 := [seq(FD[j1], j = 1 .. 2^(k-1))]; FDM1[j1] := `assuming`([simplify(combine(DiagonalMatrix(B1)), radical)], [nu > 0]) end do

Error, invalid matrix/vector product

"r:=2^(k-1)*M; print("Gegenbauer wavelets vector having order" = 2^(k-1)*M); psix :=simplify(combine(Vector[column](r, i-> psipsi[i])),radical) assuming nu>0; psit:=eval(psix,x=t); print("Unkown vector need to determine with order" = 2^(k-1)*M); C := Matrix(r,r,symbol=a); print("Trial solution of GWM when k" = k,"and M" = M); u:=simplify(combine((Transpose(psix).C.psit)),radical) assuming nu>0; for j1 from 1 by 1 while j1<=MFD do; B[j1]:=proc(M) Matrix(M, M, {seq((i,i)=((i-1)!)/(GAMMA(i- `&vartheta;`||j1)), i=j1+1..M)}); end proc: end do: for j1 from 1 by 1 while j1<=MFD do; FD[j1]:=B[j1](M) end do; for j1 from 1 by 1 while j1<=MFD do; B1:=[seq( FD[j1], j=1..2^(k-1))]: FDM1[j1]:=simplify(combine(DiagonalMatrix(B1)),radical) assuming nu>0; end do:"

  
> 

F1 := Matrix(M, M, 0); printlevel := 3; for r from 2 while r <= M do for s while s <= r-1 do if type(r+s, odd) then F1[r, s] := 2^k*sqrt(GAMMA(s-1+2*nu)*GAMMA(r)*(r-1+nu))*(2*s-2+2*nu)/sqrt(GAMMA(r-1+2*nu)*GAMMA(s)*(s-1+nu)) end if end do end do; F := F1; A := [seq(%, j = 1 .. 2^(k-1))]; D1 := simplify(DiagonalMatrix(A)); print("Operational matrix of first-order derivative"); FOD := `assuming`([simplify(combine(D1), radical)], [nu > 0]); print("Operational matrix of second-order derivative"); SOD := `assuming`([simplify(combine(D1^2), radical)], [nu > 0]); for n while n <= 2^(k-1) do for m from 0 while m <= M-1 do for j from 0 while j <= M-1 do Omega[m, j] := 2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+nu)*GAMMA(nu)^2/(Pi*2^(1-2*nu)*GAMMA(j+2*nu)))*(sum((-1)^i1*GAMMA(j-i1+nu)*2^(j-2*i1)*(sum((1/2)*binomial(m, l)*(2*n-1)^(m-l)*(1+(-1)^(j-2*i1+l))*GAMMA((1/2)*j-i1+(1/2)*l+1/2)*GAMMA(nu+1/2)/GAMMA(nu+1+(1/2)*j-i1+(1/2)*l), l = 0 .. m))/(GAMMA(nu)*factorial(i1)*factorial(j-2*i1)), i1 = 0 .. floor((1/2)*j)))/2^(k*(m+1)) end do end do; A1[n] := `assuming`([simplify(combine(Matrix(M, proc (i, j) options operator, arrow; Omega[i-1, j-1] end proc)), radical)], [nu > 0]) end do

"Operational matrix of first-order derivative"

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 4, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 8, (3, 3) = 0})

  

"Operational matrix of second-order derivative"

  

Matrix(%id = 18446744074370998870)

 (2)
> 

B := [seq(A1[j], j = 1 .. 2^(k-1))]; P := `assuming`([simplify(combine(DiagonalMatrix(B)), radical)], [nu > 0]); for j11 while j11 <= MFD do print("Operational matrix of variable-order derivative when &nu; &in;"*[j11-1, j11]); FDM[j11] := t^(-`&vartheta;` || j11)*(MatrixInverse(P).FDM1[j11].P) end do; FDM[0] := eval(FDM[1], `&vartheta;1` = `&vartheta;0`)

Error, invalid matrix/vector product

" B:=[seq( A1[j], j=1..2^(k-1))]: P:=simplify(combine(DiagonalMatrix(B)),radical) assuming nu>0: for j11 from 1 by 1 while j11<=MFD do; print("Operational matrix of variable-order derivative when &nu; in " [j11-1,j11],); FDM[j11]:=t^(-`&vartheta;`||j11)*((MatrixInverse(P).FDM1[j11]).P); end do; FDM[0]:=eval(FDM[1],`&vartheta;1`=`&vartheta;0`); "

  
> 

Transpose(FDM[0])

FDM[0]

 (3)
> 

 

Download Kitonums.mw

I think(?) these errors are associated with various "features" of 2-D math input. I generally consider debugging the latter as an unexplainable "black art". Hoiwever the following at least executes with no errors

> 

restart; with(LinearAlgebra); kernelopts(version); printlevel := 2; M := 3; k := 1; MFD := 2; `&vartheta;0` := 4/3; `&vartheta;1` := 3/4; `&vartheta;2` := 3/2; Delta := 2^(k-1)*M; nu := 1; for i while i <= 2^(k-1) do for j from 0 while j <= M-1 do `&psi;&psi;`[M*i+j-M+1] := simplify(2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+nu)*GAMMA(nu)^2/(Pi*2^(1-2*nu)*GAMMA(j+2*nu)))*(sum((-1)^i1*GAMMA(j-i1+nu)*(2*(2^k*x-2*i+1))^(j-2*i1)/(GAMMA(nu)*factorial(i1)*factorial(j-2*i1)), i1 = 0 .. floor((1/2)*j)))) end do end do; r := 2^(k-1)*M; print("Gegenbauer wavelets vector having order" = 2^(k-1)*M); `&psi;x` := `assuming`([simplify(combine(Vector[column](r, proc (i) options operator, arrow; `&psi;&psi;`[i] end proc)), radical)], [nu > 0]); `&psi;t` := eval(`&psi;x`, x = t); print("Unkown vector need to determine with order" = 2^(k-1)*M); C := Matrix(r, r, symbol = a); print("Trial solution of GWM when k" = k, "and M" = M); u := VectorMatrixMultiply(Transpose(`&psi;t`), MatrixVectorMultiply(C, `&psi;x`)); for j1 while j1 <= MFD do B[j1] := proc (M) Matrix(M, M, {seq((i, i) = factorial(i-1)/GAMMA(i-`&vartheta;` || j1), i = j1+1 .. M)}) end proc end do; for j1 while j1 <= MFD do FD[j1] := B[j1](M) end do; for j1 while j1 <= MFD do B1 := [seq(FD[j1], j = 1 .. 2^(k-1))]; FDM1[j1] := `assuming`([simplify(combine(DiagonalMatrix(B1)), radical)], [nu > 0]) end do; F1 := Matrix(M, M, 0); printlevel := 3; for r from 2 while r <= M do for s while s <= r-1 do if type(r+s, odd) then F1[r, s] := 2^k*sqrt(GAMMA(s-1+2*nu)*GAMMA(r)*(r-1+nu))*(2*s-2+2*nu)/sqrt(GAMMA(r-1+2*nu)*GAMMA(s)*(s-1+nu)) end if end do end do; F := F1; A := [seq(%, j = 1 .. 2^(k-1))]; D1 := simplify(DiagonalMatrix(A)); print("Operational matrix of first-order derivative"); FOD := `assuming`([simplify(combine(D1), radical)], [nu > 0]); print("Operational matrix of second-order derivative"); SOD := `assuming`([simplify(combine(D1^2), radical)], [nu > 0]); for n while n <= 2^(k-1) do for m from 0 while m <= M-1 do for j from 0 while j <= M-1 do Omega[m, j] := 2^((1/2)*k)*sqrt(GAMMA(j+1)*(j+nu)*GAMMA(nu)^2/(Pi*2^(1-2*nu)*GAMMA(j+2*nu)))*(sum((-1)^i1*GAMMA(j-i1+nu)*2^(j-2*i1)*(sum((1/2)*binomial(m, l)*(2*n-1)^(m-l)*(1+(-1)^(j-2*i1+l))*GAMMA((1/2)*j-i1+(1/2)*l+1/2)*GAMMA(nu+1/2)/GAMMA(nu+1+(1/2)*j-i1+(1/2)*l), l = 0 .. m))/(GAMMA(nu)*factorial(i1)*factorial(j-2*i1)), i1 = 0 .. floor((1/2)*j)))/2^(k*(m+1)) end do end do; A1[n] := `assuming`([simplify(combine(Matrix(M, proc (i, j) options operator, arrow; Omega[i-1, j-1] end proc)), radical)], [nu > 0]) end do; B := [seq(A1[j], j = 1 .. 2^(k-1))]; P := `assuming`([simplify(combine(DiagonalMatrix(B)), radical)], [nu > 0])

for j11 while j11 <= MFD do print("Operational matrix of variable-order derivative when &nu; &in;"*[j11-1, j11]); FDM[j11] := MatrixMatrixMultiply(MatrixInverse(P), MatrixMatrixMultiply(FDM1[j11], P)) end do
FDM[0] := eval(FDM[1], `&vartheta;1` = `&vartheta;0`)

`Maple 2018.1, X86 64 WINDOWS, Jun 8 2018, Build ID 1321769`

  

r := 3

  

"Gegenbauer wavelets vector having order" = 3

  

Vector(3, {(1) = 2/sqrt(Pi), (2) = (8*x-4)/sqrt(Pi), (3) = (32*x^2-32*x+6)/sqrt(Pi)})

  

Vector(3, {(1) = 2/sqrt(Pi), (2) = (8*t-4)/sqrt(Pi), (3) = (32*t^2-32*t+6)/sqrt(Pi)})

  

"Unkown vector need to determine with order" = 3

  

Matrix(3, 3, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (1, 3) = a[1, 3], (2, 1) = a[2, 1], (2, 2) = a[2, 2], (2, 3) = a[2, 3], (3, 1) = a[3, 1], (3, 2) = a[3, 2], (3, 3) = a[3, 3]})

  

"Trial solution of GWM when k" = 1, "and M" = 3

  

u := (2*(2*a[1, 1]/sqrt(Pi)+a[1, 2]*(8*x-4)/sqrt(Pi)+a[1, 3]*(32*x^2-32*x+6)/sqrt(Pi)))/sqrt(Pi)+(8*t-4)*(2*a[2, 1]/sqrt(Pi)+a[2, 2]*(8*x-4)/sqrt(Pi)+a[2, 3]*(32*x^2-32*x+6)/sqrt(Pi))/sqrt(Pi)+(32*t^2-32*t+6)*(2*a[3, 1]/sqrt(Pi)+a[3, 2]*(8*x-4)/sqrt(Pi)+a[3, 3]*(32*x^2-32*x+6)/sqrt(Pi))/sqrt(Pi)

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 2*sqrt(2)*GAMMA(3/4)/Pi, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = (16/5)*sqrt(2)*GAMMA(3/4)/Pi})

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 4/sqrt(Pi)})

  

"Operational matrix of first-order derivative"

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 4, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 8, (3, 3) = 0})

  

"Operational matrix of second-order derivative"

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 32, (3, 2) = 0, (3, 3) = 0})

  

"Operational matrix of variable-order derivative when &nu; &in;"*[0, 1]

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 4*sqrt(2)*GAMMA(3/4)/Pi, (2, 2) = 2*sqrt(2)*GAMMA(3/4)/Pi, (2, 3) = 0, (3, 1) = 0, (3, 2) = (24/5)*sqrt(2)*GAMMA(3/4)/Pi, (3, 3) = (16/5)*sqrt(2)*GAMMA(3/4)/Pi})

  

"Operational matrix of variable-order derivative when &nu; &in;"*[1, 2]

  

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 20/sqrt(Pi), (3, 2) = 16/sqrt(Pi), (3, 3) = 4/sqrt(Pi)})

  

Matrix(%id = 18446744074446164086)

 (1)
> 

Transpose(FDM[0])

Matrix(%id = 18446744074446167686)

 (2)
> 

 

Download mine.mw

The difference between the two is mainly focussed around the operators used for matrix multiplication - and I freely admit that I don't know why your original is non-functional - cos it should work!! On the other hand, I kind of assume that any use of 2-D input mode will lead to "issues" - which is why I never use it!

Well...

tomleslie 13876
July 06 2018
0 0

@Bosco Emmanuel 

I discovered something which doesn't work in "early" versions of label - the line

l:= min[index](T[..,n]):

which returns the 'index' of the minimum value. Not 100% sure when this became possible, but it doesn't work in Maple 18 \(or presumably anything earlier) - the [index] directive just isn't availabe.

I have modified my previous code to find the index of the minimum "the hard way", which is actually closer to your original code. The attached now functions in the earliest Maple release I have available - but this is still a lot later than Maple 12, so I can't guarantee anything :-(

randStuff4.mws

The multiple dots issue is apparently a "known" issue with the 2-D parser - although not known to me! I posted this as a separate question on MaplePrimes, and apparentlyMaple have known about this for ~8 years - but haven't got around to fixing it!

The mechanics of obtaining "truly random" numbers is a deep and mysterious science. Maple (like other software packages) will repeatedly generate the same set of random numbers on re-execution of the same code. Very useful for debug purposes, but not really random - so one can override this behaviour to produce a different set of random nmbers each time the same code is re-executed It is also true that if you define the same random variable, then the actual samples whihc you get will depend precisely on the way you generate th samples. There is an execution group in the attached worksheet whihc illustrates this

 

 

Which begs the obvious question...

tomleslie 13876
July 06 2018
0 0

@9009134 

What is Nxx,x?

I can see two possibilities for further progress

  1. I type these equations into Maple, having no understanding of the functional dependencies in the notation which you supply
  2. You type these equations into Maple, showing explicit functional dependencies  because you know exactly what these are

Now the first of these would take me about two minutes, but would probably be wrong because I don't understand the notation.

The second of these would take you about two minutes and (presumably) would be correct.

Which of these seems to be a better option?

consider my other comments...

tomleslie 13876
July 06 2018
0 0

@9009134 

specifically

  1. do you have six or seven dependent functions?
  2. what is meant by the subscript xx,x?

Probably not...

tomleslie 13876
July 06 2018
0 0

if I am interpreting your notation correctly

  1. You have three equations
  2. You appear to have six dependent variables, N, u__0, w__0, nu__0, M, q(x, y): maybe seven, depending on the interpretation of 'N' with the caret. As a general rule the number of equations and number of dependent variables should match,
  3. Not sure of the interpretation of subscripts with commas: I assume that Nxx is diff(N(x,y,t),x,x) and so on: but I don't understand what is meant by Nxx,x and similar structures
  4. Since three independent variables (x,y,t) are implied by your notation, you almost certainly won't be able to use numerical methods. Maple's pdsolve(pdsys, numeric) will only handle two independent variables

Where...

tomleslie 13876
July 05 2018
0 0

@vv 

is the 'spaces' requirement specified in Maple help?

As I stated...

tomleslie 13876
July 05 2018
0 0

@inteli 

It is trivial to select pretty much any combination of rows/columns from a matrix in Maple. All you need is a list which is a combination of range and individual entries. So M1B[.., [4, 1..2, 6].will select column 4, colums 1nd 2 and column 6 in that order

However you should be aware that exporting these to specific (non-contiguous) row/column locations in an Excel spreadsheet is a little more painful. Don't get me wrong - it can be done. All I/you need to know is which column of the matrix I/you want in which column of the Excel spreadsheet. (Identical argument applies to rows)

There are two ways to achieve this

  1. You provide a specification of which matrix columns go to which Excel columns - so (with the above example),  you might want
    1. M1B column 4 to Excel spreadsheet column H
    2. M1B column 1 to Excel spreadsheet column AZ
    3. M1B column 2 to Excel spreadsheet column BB
    4. M1B column 6 to Excel spreadsheet solumn D
    5. or you might want some entirely differnt mapping
    6. If you do want some more-or-less random selection of Maple matrix columns to be inserted into some more-or-less random Excel column locations, then selecting the former according to something like M1B[.., [4, 1..2, 6], is pretty much a waste of time
    7. If you do have a map of Maple matrix columns to Excel columns then supply it and someone here will do it for you
  2. Alternatively ou could try to understand the  process of selecting matrix columns and exporting these to specified Excel locations. Then you will be able to pick any of the former you want and place then anywhere in the latter you desire

Which of the above seems the better alternative?

Not your fault...

tomleslie 13876
July 05 2018
0 0

@Bosco Emmanuel 

I just had another look at this: The typo I was complaining about, ie

TC:= Array(0...1001, fill=0),

with three '.' characters in the range specification, exists in the 1-D input file which I provided!

I didn't mean to do this, because I know (honest) that the range specification should be 0..1001, with two '.' characters. The above was an error on my part or the '.' key "bounced", or something.

Interestingly: in the Maple help, it states

Note that more than two dots in succession are also parsed as the range (..) operator.

I have just run some tests, and the Maple help is half-correct:

  1. if one is using Maple's 1-D input, then two or more successive '.' characters is interpreted as the range (..) character, which is why no error/problem occurred, in this section of the provided code.
  2. if one is using Maple's 2-D input then two successive '.' characters is fine, but more than two is not (necessarily) interpreted as a range specification. It *appears* to be interpreted as a range specification if the number of '.' characters is even; but if the number of '.' characters is odd, then the final '.' character appears to be interpreted as a decimal point in front of the final number! So 1...10, *seems* to be interpreted as 1..(.10): then (.10) is interpreted as the "decimal" number 0.1. Since Array indices have to be integers, Maple 'rounds' this to 0, and the range specification becomes 1..0. At this point an Array of the "wroing" size has been defined, and causes subsequent errors.

I have fixed this issue in the attached - ie there are now no range definitions using more than two '..' characters

> 

  restart;
  with(Statistics):
#
# Initialise a few things
#
  p:= [3, 1, 0.5, 0.3, 1, 1.5]:
  RT:= Array(1..3, p[4..6]):
  TC[0]:= 0.:
  N:= Array(1..3, [0, 0, 0]):
  T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001),
                         Sample(Exponential(1/p[2]), 1001),
                         Sample(Exponential(1/p[3]), 1001)
                       ],
              scan=columns
            ):
#
# Do the calculation
#
  for n to 1001 do
      l:= min[index](T[..,n]):
      N[l]:= N[l]+1;
      B[N[l], l]:= T[l,n]+RT[l]; # changed T[1,n] to T[l,n]
      TC[n]:= TC[n-1]+B[N[l], l]
   end do:
   z:= ImportMatrix("C:/Users/TomLeslie/DeskTop/Tom.txt"):
#
# Get the result
#
   Error:= add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000);

Error := 253.811587911532

 (1)
> 

###############################################################
#
# Define procedure which does the same
# calculation as the above
#
  restart:
  lc:=proc( p::list )
          #
          # Declare/initialise local variables
          #
            local RT:= Array(1..3, p[4..6]),
                  TC:= Array(0..1001, fill=0),
                  N:= Array(1..3, [0, 0, 0]),
                  T:= Matrix( 3, 1001, [ Sample(Exponential(1/p[1]), 1001),
                                         Sample(Exponential(1/p[2]), 1001),
                                         Sample(Exponential(1/p[3]), 1001)
                                       ],
                              scan=columns
                            ),
                  z:= ImportMatrix("C:/Users/TomLeslie/DeskTop/Tom.txt"),
                  n, l, B;
          #
          # Ensure Statistics package is loaded
          #
            uses Statistics:
          #
          # Do calculations
          #
            for n to 1001 do
                l:= min[index](T[..,n]):
                N[l]:= N[l]+1;
                B[N[l], l]:= T[l, n]+RT[l]; # changed T[1,n] to T[l,n]
                TC[n]:= TC[n-1]+B[N[l], l]
            end do:
          #
          # Return result
          #
            return add((z[i,1]-(TC[i+1]-TC[i]))^2, i = 1 .. 1000);
       end proc:
#
# Run the calculation with the same data
# as used in the above
#
  Error:=lc([3, 1, 0.5, 0.3, 1, 1.5])

Error := 253.811587911532

 (2)
> 

############################################
# A 'Quirk' of sample generation
#
# Note that in the following two methods
# the 'first sample always agrees, but
# subsequent samples do not
#
  restart;
  with(Statistics):
  Sample(Exponential(1/0.3), 5);
#
  restart;
  with(Statistics):
  X:=Exponential(1/0.3):
  for i from 1 by 1 to 5 do
      Sample(X,1)[1];
  od;

Vector[row](5, {(1) = 13.726161926962995, (2) = 6.2840096650094015, (3) = 6.516070658333388, (4) = .26198922234730954, (5) = .6134828406538144})

13.7261619269630

4.78315819740088

1.51064580852117

11.1102902544401

1.52147762430717

 (3)
> 

restart; with(Statistics); p := [3, 1, .5, .3, 1, 1.5]; RT := Array(1 .. 3, p[4 .. 6]); TC[0] := 0.; N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Sample(Exponential(1/p[1]), 1001), Sample(Exponential(1/p[2]), 1001), Sample(Exponential(1/p[3]), 1001)], scan = columns); for n to 1001 do l := min[index](T[() .. (), n]); N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; z := ImportMatrix("C:/Users/TomLeslie/Desktop/Tom.txt"); Error := add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000)

Error := 253.811587911532

 (4)
> 

restart; with(Statistics); lc := proc (p::list) local RT, TC, N, T, z, n, l, B; RT := Array(1 .. 3, p[4 .. 6]); TC := Array(0 .. 1001, fill = 0); N := Array(1 .. 3, [0, 0, 0]); T := Matrix(3, 1001, [Statistics:-Sample(Exponential(1/p[1]), 1001), Statistics:-Sample(Exponential(1/p[2]), 1001), Statistics:-Sample(Exponential(1/p[3]), 1001)], scan = columns); z := ImportMatrix("C:/Users/TomLeslie/Desktop/Tom.txt"); for n to 1001 do l := min[index](T[() .. (), n]); N[l] := N[l]+1; B[N[l], l] := T[l, n]+RT[l]; TC[n] := TC[n-1]+B[N[l], l] end do; return add((z[i, 1]-TC[i+1]+TC[i])^2, i = 1 .. 1000) end proc; Error := lc([3, 1, .5, .3, 1, 1.5])

Error := 253.811587911532

 (5)
> 

 

Download randStuff3.mws

 

 

 

More an observation really...

tomleslie 13876
July 04 2018
0 0

@Mariusz Iwaniuk 

The code  aquestion_(2).mw works just fine in Maple2018.1. The final execution group produces real numbers. However the OP has stated that (s)he is using Maple 18, and in this version, the code you supply produces several "Float(undefined)" answers. See the attached for the results of running your code in Maple 18

I'm still trying to figure out why!

> 

with(LinearAlgebra):
> 

> 

r[1] := int(VectorCalculus:-`*`(VectorCalculus:-`*`(VectorCalculus:-`*`(lambda[1], alpha^2), exp(VectorCalculus:-`*`(lambda[1], Z))), 1/VectorCalculus:-`*`(VectorCalculus:-`+`(exp(VectorCalculus:-`*`(lambda[1], Z)), VectorCalculus:-`-`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(alpha))))^2, VectorCalculus:-`+`(exp(VectorCalculus:-`*`(lambda[2], Z)), VectorCalculus:-`-`(VectorCalculus:-`+`(1, VectorCalculus:-`-`(alpha)))))), Z = 0 .. infinity);

int(lambda[1]*alpha^2*exp(lambda[1]*Z)/((exp(lambda[1]*Z)-1+alpha)^2*(exp(lambda[2]*Z)-1+alpha)), Z = 0 .. infinity)

 (1)
> 

NULL

``

> 

NULL
> 

> 

R[1] := Expand(diff(r[1], lambda[1])):
> 

> 

lambda[1] := VectorCalculus:-`*`(3, 1/10):
> 

aa[1] := 0:
> 

r_r[1] := Re(simplify(r[1]))

(10/27)*3^(5/6)*arctan(3^(5/6)/(-2+3^(1/3)))+(2/9)*3^(1/6)*arctan(3^(5/6)/(-2+3^(1/3)))+(2/9)*3^(1/2)*Pi+(1/27)*3^(2/3)*ln(3^(2/3)-3^(1/3)+1)-(2/27)*3^(2/3)*ln(1+3^(1/3))-(5/27)*3^(1/3)*ln(3^(2/3)-3^(1/3)+1)+(10/27)*3^(1/3)*ln(1+3^(1/3))+2*ln(2)-(2/3)*ln(3^(2/3)-3^(1/3)+1)-(2/3)*ln(1+3^(1/3))+1/3

 (2)
> 

R_R[1] := Re(evalf(R[1]));

Float(undefined)

  

-5.223468552

  

2.445216045

  

Float(undefined)

  

Float(undefined)

  

2.445216045

  

4.417156107

  

Float(undefined)

  

0.886210077e-2

  

Float(undefined)

  

Float(undefined)

  

-0.165120849e-2

 (3)
> 

NULL
> 

 


 

Download someUndef.mw

 

First 89 90 91 92 93 94 95 Last Page 91 of 207