## 5023 Reputation

9 years, 218 days

## Duplication - apparently?...

This appears to be a duplication of the question posted here

https://www.mapleprimes.com/questions/227124-Determination-Of-Laplacian-In-A-New-Form#comment259521

under a different username - ie rather than torabi

One odd thing about that previou thread is that my final response appears to have gone missing???

In that thread I had erroneously claimed that one could not compute the Gradient() or Laplacian() of a vector function, since these required scalar arguments.

My last post (the one that has gone missing) admitted that this assertion was an error and corrected it. My worksheet accompanying that final post shows that (using the VectorCalculus() package) one can compute the Laplacian() of a vector function, (the result is a vector), but one cannot compute the Gradient() of a vector function (technically the result is a second order tensor, aka a matrix). I pointed out that (with care) one *ought* be able to use the transpose of the Jacobian as a proxy for the Gradient() of a vector function.

This "missing" worksheet is reposted here

 > restart;   with(VectorCalculus): # # Define coordinates #   X := a*sinh(tau)*cos(phi)/(cosh(tau) - cos(sigma)):   Y := a*sinh(tau)*sin(phi)/(cosh(tau) - cos(sigma)):   Z := a*sin(sigma)/(cosh(tau) - cos(sigma)): # # Set up coordinate system. NB this system which the # OP calls "toroidal" is slightly different from the # Maple's in-built toroidal system - hence the # coordinate system name #   AddCoordinates( OPToroid[ sigma,                             tau,                             phi                           ],                  [ X, Y, Z ]): #
 > SetCoordinates(OPToroid[sigma, tau, phi]): V:=VectorField( );
 (1)
 > simplify(expand(Laplacian(V)));
 (2)
 > LinearAlgebra:-Transpose(Jacobian(V));
 (3)
 >

## I wouldn't consider...

either of the animations in the attached as "consuming a lot of resources"

 > restart;   with(plots):   P:=[1,1]:   Q:=[2,2]:   pltOpts:=symbol=solidcircle, symbolsize=20, color=red:   nf:=50:   animate( pointplot,            [ (1-i/nf)*~P+(i/nf)*~Q,              pltOpts,              title="Point Only"            ],            i=0..nf,            frames=nf          );   animate( pointplot,            [ [ 'seq                 ( (1-j/nf)*~P+(j/nf)*~Q,                   j=0..i                 )'              ],              pltOpts,              title="Trail"            ],            i=0..nf,            frames=nf          );
 >
 >

## More of an observation really...

No problem executing any of your code in Maple 2019 - see the attached.

The output of your command sin(x) is interesting. Since 'x' is undefined, Maple ought to return 'sin(x)'. The only occasions reported here where random numeric values are returned in such a situation, is when an illegal (aka "cracked") version of Maple is being used.

Assuming you haev a "paid for" version of Maple with a valid activation code, I can only suggest that you contact Maple support

 > plot(sin(x), x = -2*Pi .. 2*Pi);
 > sin(x);
 (1)
 > plot3d(x*exp(-x^2 - y^2), x = -2 .. 2, y = -2 .. 2, color = x);
 > with(GraphTheory): G := Graph({{a, b}, {a, c}, {b, c}}); DrawGraph(G)
 >

## Something strange?...

1. Carl is correct about the name conflict (fixed in the attached)
2. Even with this "fix", the numeric solver fails on the first ODE system, because the "initial Newton iteration is not converging". There are several possible causes of this problem, which may be worth investigating, but
3. The information you require ( ie y(t) ), can be obtained by solving only the first two ODEs in the system with their associated boundary conditions. The other five ODEs in the system (togther with their associated boundary conditions) appear to be "superfluous". This strikes me as odd - why would you define them, if you don't need them??
4. It is possible that the non-convergence problem and the "superfluous" equations are related in some way, but I think you should carefully consider why your system contains equations you don't seem to need. I'm geussing that either I have got something wrong somewhere - or you have!
5. So in the attached I have solved a "reduced" version of your first ODE system. It *may* be possible to use the solutions for x(t) and y(t) to solve the other equations in this system (if indeed it is necessary)

Anyway, for what it is worth, the attached shows the solutions

 > restart: with(plots): r := 3: r__1 := 3: k := 10: a := 0.2e-1: b := 0.1e-1:  c := 0.1e-1: beta := 0.3e-1: alpha := 0.3e-1: m := 0.5e-1: z := 40: q := 5: p := 100: T := 3: sigma := 0.1e-1: k__1 := 10: rho := 0.5e-1:  u[1] := min(max(0, z), 1):  z := (a*m*k*lambda[2](t)*x(t)*y(t)-lambda[1](t)*r*(1+b*x(t)+c*y(t))*x(t)*x(t))/(z*k*(1+b*x(t)+c*y(t))):  u[2] := min(max(0, q), 1):  q := -lambda[1](t)*beta*x(t)*s(t)/q:  u[3] := min(max(0, p), 1):  p := -(r__1*lambda[3](t)*s(t)*s(t))/(p*k__1):  sys := diff(x(t), t) = r*x(t)*(1-(1-u[1])*x(t)/k)-a*m*x(t)*y(t)/(1+b*x(t)+c*y(t))-beta*(1-u[2])*x(t)*s(t),         diff(y(t), t) = -alpha*y(t)+a*m*x(t)*y(t)/(1+b*x(t)+c*y(t)),         diff(s(t), t) = sigma*s(t)+r__1*s(t)*(1-(1-u[3])*s(t)/k__1)-rho*s(t)*y(t),         diff(lambda[1](t), t) = -lambda[1](t)*(r-2*r*(1-u[1])*x(t)/k-a*y(t)*(1+c*y(t))/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t)))-beta*(1-u[2])*s(t))-lambda[2](t)*a*m*(1-u[1])*(1+c*y(t))*y(t)/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t))),         diff(lambda[2](t), t) = -lambda[1](t)*a*x(t)*(1+b*x(t))/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t)))+lambda[2](t)*(-alpha+(a*m*(1-u[1]) . (1+b*x(t)))*x(t)/((1+b*x(t)+c*y(t))*(1+b*x(t)+c*y(t))))+lambda[3](t)*rho*s(t),        diff(lambda[3](t), t) = lambda[1](t)*beta*(1-u[2])*x(t)-lambda[1](t)*(r__1-2*r__1*(1-u[3])*s(t)/k__1-sigma-rho*y(t)),        x(0) = 10, y(0) = 20, s(0) = 10, lambda[1](T) = 0, lambda[2](T) = 0, lambda[3](T) = 0: #p1 := dsolve({sys}, type = numeric):   p1:=dsolve([sys[1..2], sys[7..8]], numeric):   p2o := odeplot( p1,                   [t, y(t)],                   0 .. 2,                   labels = ["Time (months)", " "*`badbiomass""spatina"""`],                   labeldirections =   [horizontal, vertical],                   color = red,                   axes = boxed                 );
 > r := 3: r__1 := 3: k := 10: a := 0.2e-1: b := 0.1e-1:   c := 0.1e-1: beta := 0.3e-1: alpha := 0.3e-1: m := 0.5e-1:   sigma := 0.1e-1: k__1 := 10: rho := 0.5e-1:   z := 40: q := 5: p := 100: T := 3:   fun := diff(x(t), t) = r*x(t)*(1-x(t)/k)-a*m*x(t)*y(t)/(1+b*x(t)+c*y(t))-beta*x(t)*s(t),          diff(y(t), t) = -alpha*y(t)+a*m*x(t)*y(t)/(1+b*x(t)+c*y(t)),          diff(s(t), t) = sigma*s(t)+r__1*s(t)*(1-s(t)/k__1)-rho*s(t)*y(t),          x(0) = 10, y(0) = 20, s(0) = 10:   p2 := dsolve({fun}, type = numeric);   p2i := odeplot( p2,                   [t, y(t)],                   0 .. 2,                   labels = ["Time(months)", " bad biomass"],                   labeldirections = [horizontal, vertical],                   axes = boxed,                   color = blue                 );   plots[display](p2i, p2o);
 >

## You might...

want to consider the patmatch() function

## Result not wrong!...

You might want to check the explanation in the example worksheet returned by

?Examples/functionaloperators

In particular the distinction between writing 1(x) in 1-D input and 1(x) in 2-D math input - because the first will return '1', and the second will return 'x'

## Search space is small...

so exhaustive enumeration is feasible. NB this approach would "scale" very badly. There are (almost certainly!) more elegant approaches which would scle better

In the attached I have made a number of assumptions, since your precise conditions are unclear/ambiguous. These are

1. A>0, B>0, C>0. Allowing any of these to be zero just means changing the start index on the "loops" in the attached
2. in the requirement that `or`(A,B,C)=2^n, the value n=0 is allowed. Maybe you wnat the minimum value of 'n' to be 1? Again this would be trivial to fix
3. when you say " A or B can be prime or 2^n  * prime", this is ambiguous. I have interpreted this condition as

or(   or(A is prime,   A is 2^n)
or(B is prime,   B is 2^n )
)

Since the problem is symmetrical, any permutation of the results obtained would also be a solution

With the above assumptions, check the attached

 > restart;
 > # # generate list of relevant primes #   primList:= [ seq                ( `if`                  ( isprime(j),                    j,                    NULL                  ),                  j=1..60                )              ]: # # generate list of relevant powers of 2 # NB 2^0 is allowed #   powList:= [ seq               ( 2^j,                 j=0..5               )             ]: # # generate list of of relevant products #   prodList:= select              ( i-> i<60,                { seq                  ( (j*~powList)[],                    j in primList                  )                }              ): # # Generate some relevant logical tests #   cond1:= z-> member(z, primList):   cond2:= z-> member(z, powList):   cond3:= z-> `or`( cond1(z), cond2(z)):
 > ########################################## # Check all combinations for Question 1 # # Initialise solution #   sol1:= [a=0, b=0, a+b=0]:   for i from 1 to 60 do       for j from 1 to 60-i do           if   `and`                ( `or`( cond1(i-j),                        cond2(i-j)                      ),                  `or`( cond3(i),                        cond3(j)                      ),                  `or`( cond1(i),                        cond1(j)                      )                )           then if   i+j>rhs(sol1[3])                then sol1:=[a=i,b=j,a+b=i+j];                fi:           fi;        od;   od;   sol1;
 (1)
 > ############################################ # Initialize solution #   sol2:= [a=0, b=0, c=0, a+b+c=0]: # # Check all combinations for Question 2 #   for i from 1 to 60 do       for j from 1 to 60-i do           for k from 1 to 60-j do               if   `and`                    ( `or`( cond1(i-j),                            cond2(i-j)                          ),                      `or`( cond1(j-k),                            cond2(j-k)                          ),                      `or`( cond1(i),                            cond1(j),                            cond1(k)                          ),                      `or`( cond3(i),                            cond3(j),                            cond3(k)                          )                    )               then if   i+j+k > rhs(sol2[4])                    then sol2:= [a=i,b=j,c=k, a+b+c=i+j];                    fi:               fi;           od;       od;   od;   sol2;
 (2)
 >

## Aaaah -Eratosthenes sieve...

This is a classic - so I dug ot some code I wrote many years ago, which is a (reasonably) efficient implementation of the Sieve of Eratosthenes (check Wikipedia for an explanation)

 > sieve:= proc( N::integer ) :: list;                 local j, k, prArr:                                 description "The sieve of Eratosthenes":               #               # Initialise the array with all even indices               # >=4 set to zero: other entries set to 1               #                 prArr:= Array                         ( 2..N,                           [ 1,                             1,                             seq                             ( '0,1',                               j = 2..N/2                             )                           ]                         ):               #               # cycle through all odd multiples (other than 1)               # of all odd indices, setting the entry to 0               #                 prArr[ [ seq                                      ( seq                            ( k*j,                              j = k..N/k, 2                            ),                            k = 3..trunc(sqrt(N)), 2                          )                        ]                      ]:= 0:               #               # Return the indices which have non-zero entries               #                 return [ 2,                          seq                          ( `if`                            ( prArr[j]=0,                              NULL,                              j                            ),                            j = 3..N, 2                          )                        ]:                     end proc:
 > sieve(10000);
 (1)
 >

## For the terms you have shown...

use the genfunc() package, as in the attached

 > restart;   with(genfunc):   p:=[2,3,4,6,8,9,10,12,14,15,16,18]; # # Produce recurrence relation #   rgf_findrecur(6, p, f, j);
 (1)
 > # # Check the above recurrence relation by trying # to regenerate the sequence #   g:=k-> if   k=1          then 2          elif k=2          then 3          elif k=3          then 4          elif k=4          then 6          else 2*g(k-1)-2*g(k-2)+2*g(k-3)-g(k-4)          fi;   [ seq     ( g(i),       i=1..numelems(p)     )   ];
 (2)
 >

## Exact numbers -maybe!...

I'm not sure the attached is 100% valid, but it does get the *same* answer for the maximal distance, which appears to be

sqrt(3)*(sqrt(2)+4)

which is obtained with  a=1, b=2/3, c=1/3 which are close(?) to the values obtained numerically.

Anyhow, FWIW, check out the attached

 > restart;   with(geom3d):   point(A, [2,0,0]):   point(B, [0,3,0]):   point(C, [0,0,6]):   assume(a<>0):   interface(showassumed=0):   v:= [a,b,c]:   point(DD, [1,1,1]):   line(l1, [DD, v]):   d1:=distance(A, l1)+distance(B, l1)+distance(C,l1); # # Just to make things easier for the maximize() function # use the squares of the individual distances #   d2:= distance(A, l1)^2+distance(B, l1)^2+distance(C,l1)^2: # # and simplify a bit #   d2:= simplify(expand(d2)); # # Get a "symbolic" maximum #   sol:= maximize(d2, a=0..1, b=0..1, c=0..1, location); # # Evaluate the "actual" distance function with these values #   simplify( eval(d1, sol[2,1,1] ), symbolic); # # Get a floating point value for the above, for comparison # with a numerical optimization #   evalf(%);   Optimization:-Maximize(d1, a=0..1, b=0..1, c=0..1);
 (1)
 >

## As far as I can tell...

you have two problems, one syntactic and one logical

Syntatic

In Maple you can run into all sorts of weird problems if you use the same name for indexed and unindexed variables - so writing anything like f(x):=sum(f[i](x),i=1..N) is just asking for trouble. Where you have done this I have capitalised the left-hand-side to produce distict names. (Obviously, I have appropriately corrected all subsequent name references in the code)

Logical

Becuase your equations are coupled, you cannot solve for f[i] i=0..N, then theta[i] i=0..N and so on sequentially. However for any value of 'i' you can simultaneously solve for the [ f[i], theta[i], u[i], w[i] ], since these only depend on previously obtained values, ie [ f[i-1], theta[i-1], u[i-1], w[i-1] ].

The attached now executes withut error, and appears to return all required quantities Be aware that some of these expressions are a bit "lengthy"

 > restart;   PDEtools[declare](f(x),theta(x),u(x),w(x), prime=x):   N:= 4;
 (1)
 > F:=      x -> local i;                 add                 ( p^i*f[i](x),                   i=0..N                 );   Theta:=  x -> local i;                 add                 ( p^i*theta[i](x),                   i = 0..N                 );   U:=      x -> local i;                 add                 ( p^i*u[i](x),                   i = 0..N                 );   W:=      x -> local i;                 add                 ( p^i*w[i](x),                   i = 0 .. N                 );
 (2)
 > HPMEq := collect            ( (1-p)*diff(F(x), x\$2)+p*(diff(F(x), x\$2)-k1*diff(F(x), x)-k2*F(x)),               p            );   HPMEr := collect            ( (1-p)*diff(Theta(x),x\$2)+p*(diff(Theta(x),x\$2)-k11*diff(Theta(x),x)+k12*diff(U(x),x)^2+k13*diff(W(x),x)^2+k14*Theta(x)),               p            );   HPMEs := collect            ( (1-p)*diff(U(x),x\$2)+p*(diff(U(x),x\$2)-R*diff(U(x),x)-A-k8*W(x)-k7*U(x)+k5*Theta(x)+k6*F(x)),               p            );   HPMEt := collect            ( (1-p)*diff(W(x),x\$2)+p*(diff(W(x),x\$2)-R*diff(W(x),x)+k9*U(x)-k10*W(x)),               p            );
 (3)
 > # # renamed equ[1] to equ1 #   for i from 0 to N do       equ1[i] := coeff(HPMEq, p, i) = 0;   end do; # # renamed equa[1] to equ2 #   for i from 0 to N do       equ2[i] := coeff(HPMEr, p, i) = 0;   end do; # # renamed equat[1] to equ3 #   for i from 0 to N do       equ3[i] := coeff(HPMEs, p, i) = 0;   end do; # # renamed equati[1] to equ4 #   for i from 0 to N do       equ4[i] := coeff(HPMEt, p, i) = 0;   end do
 (4)
 > # # Sequentially solve the system(s) #   for i from 0 to N do       cons:= f[i](-1) = 1, f[i](1) = 1,              theta[i](-1)=1, theta[i](1)=1,              u[i](-1)=1, u[i](1)=1,              w[i](-1)=1, w[i](1)=1;       assign(dsolve( [ equ1[i], equ2[i], equ3[i], equ4[i], cons ]));   od:
 > # # What happens next?? # # OP can access any of the "individual" functions, # for example theta[2](x) or w[3](x), just by entering # these expressions, as in #   collect(theta[2](x), x);   collect(w[3](x), x); # # Or one can get the "sum" terms for these functions # just by entering F(x), Theta(x), U(x) or W(x). These # expresssions are a bit lengthy # # For example, U(x) is shown below #   collect(U(x), x);
 (5)
 >

## Problem with using ->...

Defining a function with the '->' notation only really allows very simple functions. In particular 'if' statements are only allowed to be of the form

if(condition, doThis, otherwiseDoThis)

'elif' clauses are not permitted. Thus your definiton of the function 'q1' is invalid. Interestingly, you have perfectly valid definition of 'q1' as a proc(), which is commented out! If I comment out the 'q1->....' statement and uncomment the qi:=proc(....) statement, then all I get are a few warnings in the calling procedure about undeclared local variables. This is easily fixed by adding these to the 'local' statement.

The result executes with no errors - obviously I have not checked whether the resulting procedure provide the answers you want, bu the attached executes with no errors.

## as in...

the attached, maybe(?)

 > restart;   interface(version);
 (1)
 > with(LinearAlgebra):   interface(rtablesize=20):   with(plots):   k:=2:M:=2:   Tm:=(t,m)-> 2*t*Tm(t,m-1)-Tm(t,m-2);   Tm(t,0):=1:   Tm(t,1):=t:   Tm2:=(tt,m)->subs(t=tt,Tm(t,m)):   alpha:=(m)->piecewise(m=0,1/sqrt(Pi),sqrt(2)/sqrt(Pi));   psii:=(n,m,x)->piecewise((n-1)/(2^(k-1)) <= x and x <= n/(2^(k-1)),   alpha(m)*(2^(k/2))*Tm2(2^k*x-2*n+1,m), 0);   psi:=(t)-> local i, j;              Array([seq(seq(psii(i,j,t),j=0..M-1),i=1...2^(k-1))] ):   for i from 1 to ((2^(k-1))*M) do       r[i]:=evalf(psi((2*i-1)/((2^k)*M))):   end do:   m:=M*(2^(k-1)):   xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));   PB_n:= Matrix( m,                  m,                  (i, j)-> `if`( i=j,                                 1,                                 `if`( j>i,                                       xi(j-1,n),                                       0                                     )                                )                 );   PB_n:=1/m^n/factorial(n+1)*PB_n;
 (2)
 >

## Actually...

the following version is probably better

 > restart; eq:=z=(r-4)/(r-2)/(r-6)+15: sol:=solve(eq, r): plot3d( sol[2], theta=0..2*Pi, z=0..40, coords=cylindrical, grid=[100,100], style=surface, axes=none);
 >