Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

About once a year Advent of Code give you a problem that is a gift if you are using a computer algebra system. This year that day was Day 13. The day 13 problem was one of crazy claw machines.  Each machine has two buttons that can be pressed to move the claw a given number of X and Y positions and a prize at a given position. A buttons cost 3 tokens to press, and B buttons cost 1 token to press, and we are asked to find the minimum number of tokens needed, IF it is possible to reach the prize. The data presented like so:

Button A: X+94, Y+34
Button B: X+22, Y+67
Prize: X=8400, Y=5400

Button A: X+26, Y+66
Button B: X+67, Y+21
Prize: X=12748, Y=12176

...

Some times the input is harder to parse than the problem is so solve, and this might be one of those cases.  I tend to reach for StringTools to take the input apart, but today the right tool is the old school C-style sscanf (after using StringSplit to split at every double linebreak).

machinesL := StringTools:-StringSplit(Trim(input), "\n\n");
machines := map(m->sscanf(m,"Button A: X+%d, Y+%d\nButton B: X+%d, Y+%d\nPrize: X=%d, Y=%d"), machinesL):

Now we have a list of claw machine parameters in the form [A_x, A_y, B_x, B_y, P_x, P_y] and we need to turn those into equations that we can solve. We want the number of A presses a, and B presses b to get the claw to the P_x, P_y position of the claw, it is simple to just write them down:

for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
end do;

Now because of the discrete nature of this problem, we need our variables a and b to be non-negative integers.  When solving this, I first reached for isolve like this:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol <> NULL then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Now, sometimes Advent of Code inputs contain a lot of hidden structure.  I wrote the code above, it worked on the sample input, so I tried it immediately on my real input (about 300 claw machines like the above) and IT WORKED.  But, you might notice that this code does not deal with a couple cases that could have appeared.  In particular, it doesn't check that the solutions are positive.  It also doesn't handle cases where there is more than one possible solution.  The former is easy to check

if sol <> NULL and eval(a,sol) >= 0 and eval(b,sol) >= 0 then

Unfortunately isolve does not handle inequalities, but you could try with solve, but it doesn't save us any checking, because we'd still have to check if the solutions are integers, so we might as well have just solved the equation and then checked if it were a nonnegative integer.

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := solve(eqn);
   if type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;
ans1 = tokens;

In the multiple solution case we get something like {a = 3 - 2*b, 0 <= b, b <= 3/2} which has some great information in it but might be hard to handle programmatically, so let's see what isolve does with those cases to see if it's easier to deal with

> eqn := { 17*a + 84*b = 7870, 34*a + 168*b = 15740 }:
> constr := { a >= 0, b>= 0 }:
> sol := isolve(eqn);
               sol := {a = 458 - 84 _Z1, b = 1 + 17 _Z1}

> constr := eval({ a >= 0, b>= 0 }, sol):
            const := {0 <= 1 + 17 _Z1, 0 <= 458 - 84 _Z1}

> obj := eval(3*a+b, sol);
                         obj := 1375 - 235 _Z1

You can see it's easy to tell if these show up in your input, since your "token" total will have the _Zn variables in it.  Now, since everything is simple and linear here, it seems like you could use solve to find the rational value of _Z1 that makes obj=0 and then take the closest integer but it's not so simple, we actually have to deal with the contraints that a and b be positive too. So, it really just makes sense to bring out the big hammer of Optimization:-Minimize which allows us to directly optimize over just the integers.  So a full solution looks like this:

tokens := 0:
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   sol := isolve(eqn);
   if sol = NULL then
      next;
   end if;
   constr := eval({ a >= 0, b>= 0 }, sol);
   obj := eval(3*a+b, sol);
   if not type(obj, constant) then
      tokens := tokens + Optimization:-Minimize( obj, constr, assume=integer )[1];
   elif andmap(evalb, constr) then
      tokens := tokens + obj; 
   end if;
end do;

But since we're bringing out the big hammer, why not just use Optimization in the first place.  The main reason is that Minimize doesn't simply return NULL when it doesn't work, instead it throws an exception, so we need to find all the exceptions that can occur and handle then with a try-catch, thus:

tokens := 0;
for m in machines do
   eqn := ({m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6]});
   try 
       sol := Optimization:-Minimize(3*a+b, eqn, 'assume'='nonnegint')[1];
       tokens := tokens + sol;
   catch "no feasible":
   end try; 
end do;
tokens;

(you can in fact omit the string in the catch: statement, but I can tell you from long experience that that is an excellent way to make your code much much harder to debug)

Alright, so how did people not using Maple solve this problem?  The easiest way to solve it, and the one used by all the cheaters scraping the website and using LLM-based code generators that auto-submit solutions to get into the Top 100, was to just check all possible a, b values in 0..100 and take the values than minimize 3*a+b when reaching the prize coordinates.  That's only feasible because the problem states the 100 is an upperbound for a and b, but it's also very fast (about 1/10 second in Maple):

tokens := 0:
for m in machines do;
sol := infinity;
for i from 0 to 100 do for j from 0 to 100 do
    if i*m[1]+j*m[3]=m[5] and i*m[2]+j*m[4]=m[6] and 3*i+j < sol then
        sol := 3*i+j;
    end if;
end do; end do;
tokens := tokens + ifelse(sol=infinity,0,sol);
end do;

It does not scale at all to part 2 (which modified everything to be bigger by about 10 trillion), and it seems that foiled all the LLM solvers. So, what solutions scaled in languages without integer equation solvers?  Well, the easiest solution is just to solve the general equation using paper and pencil


And you can just hard code that formula in, check that it gives integer values and compute the tokens. As long as you get unique solutions, that looks something like this

solveit := proc(m)
local asol := m[4]*m[5] - m[3]*m[6];
local bsol := m[1]*m[6] - m[2]*m[5];
local deno := m[1]*m[4] - m[2]*m[3];
if deno = 0 then return -2^63; end if; # multiple solution case - not handled
if asol mod deno = 0 and bsol mod deno = 0
   and (   ( deno>=0 and asol>=0 and bsol>=0 ) 
        or ( deno<=0 and asol<=0 and bsol<=0 ) )
then

    return 3*asol/deno+bsol/deno;
else
    return 0;
end if;
end proc:

Which if you have this in Maple, you can impress your friends by auto generating solutions in other languages. Here, for your FORTRAN friends

> CodeGeneration:-Fortran(solveit);

Warning, the following variable name replacements were made: solveit -> cg
       integer function cg (m)
        doubleprecision m(*)
        integer asol
        integer bsol
        integer deno
        asol = int(-m(3) * m(6) + m(4) * m(5))
        bsol = int(m(1) * m(6) - m(2) * m(5))
        deno = int(m(1) * m(4) - m(2) * m(3))
        if (deno .eq. 0) then
          cg = -9223372036854775808
          return
        end if
        if (mod(asol, deno) .eq. 0 .and. mod(bsol, deno) .eq. 0 .and. (0
     # .le. deno .and. 0 .le. asol .and. 0 .le. bsol .or. deno .le. 0 .a
     #nd. asol .le. 0 .and. bsol .le. 0)) then
          cg = 3 * asol / deno + bsol / deno
          return
        else
          cg = 0
          return
        end if
      end

Another way that you might solve this without solve is to use a linear algebra library to solve the linear system.  It works even if you only have a numeric solver, but you have to be careful about checking for integers:

tokens := 0:
for m in machines do
    sol := LinearAlgebra:-LinearSolve(
               Matrix(1..2,1..2,[m[[1,3]],m[[2,4]]], datatype=float), 
               Vector(m[5..6], datatype=float));
    if abs(sol[1]-round(sol[1])) < 10^(-8) and abs(sol[2]-round(sol[2])) < 10^(-8)
       and sol[1] >= 0 and sol[2] >= 0
    then
       tokens := tokens + 3*sol[1]+sol[2];
    end if;
end do;

Finally, a lot of people solved this sort of thing with the Z3 Theorem prover from Microsoft research which is also way more than you need, but it mostly just uses SMTLIB, which we also have in a library for in Maple, and it can just be used in place of solve

tokens := 0;
for m in machines do
   eqn := {m[1]*a+m[3]*b=m[5], m[2]*a+m[4]*b=m[6], a>=0, b>=0};
   sol := SMTLIB:-Satisfy(eqn) assuming a::nonnegint, b::nonnegint;
   if sol <> NULL and type(eval(a,sol), integer) and type(eval(b,sol),integer) then
      tokens := tokens + eval(3*a+b, sol);
   end if;
end do;

Notice that Satisfy handled the multiple solution case just by choosing one of the many solutions. It is possible to get SMTLib to optimize but it is slightly more involved, and this post is already too long. This time, I've put all this work in worksheet: Day13-Primes.mw

Hi!

So I like to check that my manual integrations and/or maple integrations are equal with each other. I normally do this using the Test Relation function.

I was working on a problem and noticed that Maple didn't evaluate the integrals being the same, even though they presumedly are.

Could anyone shed some light on why I get this inequality?

Thanks in advance!

mapleintvsmanualint.mw

restartNULL

dn/dt = -r__S*V

 

We can define n as C*V 

dC*V/dt = -r__S*V

 

We can define the concentration C as S  

dS/dt = -r__S

 

where -r__S = V__max*[S]/(1+K__1*[S]+K__2*[S^2])

  

ds/dt = -V__max*[S]/(1+K__1*[S]+K__2*[S^2])=

 

1/dt = -V__max*[S]/((1+K__1*[S]+K__2*[S^2])*ds)

 

dt = (1+k__1*[S]+K__2*[S^2])*ds/(V__max*[S])

 

"&DifferentialD;t = 1/(`V__max`*[S])+`k__1`/(`V__max`)+(`K__2`*[S])/(`V__max`)*&DifferentialD;s"

 

int(1/(V__max*S)+K__1/V__max+K__2*S/V__max, S = S .. S__0)

 

`assuming`([simplify(combine*(int(1/(V__max*S)+K__1/V__max+K__2*S/V__max, S = S__ .. S__0)), size)], [S > 0, S__0 > S__])

combine*piecewise(And(0 < S__0, S__ < 0), undefined, (1/2)*(K__2*S__0^2-K__2*S__^2+2*S__0*K__1-2*K__1*S__+2*ln(S__0)-2*ln(S__))/V__max)

(1)

 

`assuming`([simplify(int(1/(V__max*S)+K__1/V__max+K__2*S/V__max, S = S__ .. S__0), size)], [S > 0, S__0 > S__])
  piecewise(And(0 < S__0, S__ < 0), undefined, (1/2)*(K__2*S__0^2-K__2*S__^2+2*S__0*K__1-2*K__1*S__+2*ln(S__0)-2*ln(S__))/V__max)NULL

 

 

maple*equation = manual*equation  NULL

(S__0^2*K__2-K__2*S^2+2*S__0*K__1-2*K__1*S+2*ln(S__0)-2*ln(S))/(2*V__max) = (ln(S__0/S)+K__1*(S__0-S)+(1/2)*(-S^2+S__0^2)*K__2)/V__max"(->)"false

   

eq1 := (S__0^2*K__2-K__2*S^2+2*S__0*K__1-2*K__1*S+2*ln(S__0)-2*ln(S))/(2*V__max)

 

eq2 := (ln(S__0/S)+K__1*(S__0-S)+(1/2)*(-S^2+S__0^2)*K__2)/V__max

 

eq1-eq2 = 0"(->)"false

 

(ln(S__0/S)+K__1*(S__0-S)+(1/2)*(-S^2+S__0^2)*K__2)/V__max = (ln(S__0/S)+K__1*(S__0-S)+(1/2)*(-S^2+S__0^2)*K__2)/V__max

 

Download mapleintvsmanualint.mw

The uploaded worksheet defines a surface.

I would like to code a ball rolling across this surface (and others) starting from an initial position on the surface and an initial velocity tangent to the surface, but I don't know how to do this.

What is the combination of physics (including gravity) and math that accomplishes this task?

Ball-rolling-on-a-surface.mw

Is there a trick to make Maple give same result below when using eval and limit?  

Attached worksheet. This comes in context of solving ode  using Laplace. Initial conditions are at zero. And need to solve for the constant of integration. 

It works when using eval, since Dirac(t) becomes Dirac(0), but when using Limit, Dirac(t) becomes zero and the _C1 is lost. I was wondering if limit should also return Dirac(0) like eval?

interface(version);

`Standard Worksheet Interface, Maple 2024.2, Windows 10, October 29 2024 Build ID 1872373`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1838 and is the same as the version installed in this computer, created 2024, December 2, 10:11 hours Pacific Time.`

restart;

e:=1/2*t+_C1*Dirac(t);

(1/2)*t+_C1*Dirac(t)

eval(e,t=0)

_C1*Dirac(0)

limit(e,t=0)

0

 

 

Download dirac_limit_dec_13_2024.mw

Some of the calculations mentioned here can be done in alternative programming languages, such as Python, C, and so on. However, I would like to reproduce exactly these graphs using Maple (without the need for programming commands, such as "if", "while", among others).

In the work I am trying to reproduce, we have "The evaluation of the influence of the inclusion of the broadband behavior of grounding systems in EMT-type programs in the evaluation of transients resulting from direct lightning strikes on transmission lines. The behavior of the grounding frequency is determined using an accurate electromagnetic model and included in the EMTP/ATP by means of an equivalent circuit derived from the Vector Fitting technique. In addition, the impact of the frequency dependence of soil parameters on the lightning performance of transmission lines is addressed." This may seem somewhat disconnected from reality for many, since it is a problem involving electrical engineering optimization.

Could someone help me reproduce these calculations? I have made little significant progress.

If you want to access the reference accounts, I'll send you the PDF

schroeder2017 [link to copyrighted material replaced by moderator]

Hello Dear Maple Users and Experts

I am running this code for N=3, but fsoolve can not work and did not give me any output. Could you help me how can I get the result?

Actually, I got the result for N=2. Exact solution is a[0]=0, a[1]=1, a[2]=1, b[1]^2=sqrt(2) and b[1]^2+b[2]^2=sqrt(3). But, for N=3, I can not receive any results from fsolve

Here is my code

restart;
Digits := 20;
L := 1;
N := 3;
alpha := 1;
xexact := t -> t^sqrt(2) + t^sqrt(3);
f := simplify(fracdiff(t^sqrt(2), t, alpha)) + simplify(fracdiff(t^sqrt(3), t, alpha));
f := unapply(f, t);
xapp := a[0] + sum(a[j]*t^sum(b[i]^2, i = 1 .. j), j = 1 .. N);
xapp := unapply(xapp, t);
xfrac := sum(a[jj]*simplify(GAMMA(sum(b[ii]^2, ii = 1 .. jj) + 1)/GAMMA(sum(b[ii]^2, ii = 1 .. jj) + 1 - alpha))*t^(sum(b[ii]^2, ii = 1 .. jj) - alpha), jj = 1 .. N);
xfrac := unapply(xfrac, t);
xfrac1 := sum(a[jj]*simplify(sum(b[ii]^2, ii = 1 .. jj)^(alpha + 1)/(sum(b[ii]^2, ii = 1 .. jj) - alpha))*t^(sum(b[ii]^2, ii = 1 .. jj) - alpha), jj = 1 .. N);
xfrac1 := unapply(xfrac1, t);
S1 := {seq(evalf(xfrac(k/(2*N)*L)) - evalf(f(k/(2*N)*L)) = 0, k = 1 .. 2*N)};
S2 := {xapp(0) = 0};
S := S1 union S2;
sol := fsolve(S);

Hi,
How can I simplify this relation(See uploaded .mw file)?
For example, the second term is simplified as: 

deltae*(1-phi0/(kappa-3/2))^(-kappa+1/2)+(1/2)*deltab*(1-sqrt(2)*sqrt(1/(m*ub^2))*sqrt(-phi0));

di1.mw

In Latest Maple 2024.2, I found that when doing z:=%  where % is result on integration, causes internal error 

          Error, unexpected result from Typesetting

But when the interface is set to standard, no such error.

This not only happen in worksheet, but also when code is run in command line!

Worksheet below.

interface(version);

`Standard Worksheet Interface, Maple 2024.2, Windows 10, October 29 2024 Build ID 1872373`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1837 and is the same as the version installed in this computer, created 2024, December 2, 10:11 hours Pacific Time.`

Example using extended

 

restart;

interface(typesetting=extended):

int(exp(-int(b(t),t))*t^4*csc(t)^2,t);

int(exp(-(int(b(t), t)))*t^4*csc(t)^2, t)

z:=%;

Error, (in Risch:-Norman) too many levels of recursion

` `

Error, unexpected result from Typesetting

 

Example using standard

 

restart;

interface(typesetting=standard):

int(exp(-int(b(t),t))*t^4*csc(t)^2,t);

int(exp(-(int(b(t), t)))*t^4*csc(t)^2, t)

z:=%;

int(exp(-(int(b(t), t)))*t^4*csc(t)^2, t)

 

 

Example using direct assignment also

 

restart;

interface(typesetting=extended):

z:=int(exp(-int(b(t),t))*t^4*csc(t)^2,t);

Error, (in Risch:-Norman) too many levels of recursion

` `

Error, unexpected result from Typesetting

Download extended_interface_causes_internal_bug_dec_13_2024.mw

ps. also reported to Maple support

I was rejected because the editor said my equation is too long. My question is: Is there a way to rewrite the equation in a more concise form? Additionally, is there a package in Maple that allows for automatic simplification or collection of terms without using specific commands? Any suggestions for addressing this issue would be appreciated.

restart

``

eq3 := -6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^3*a[4]+(6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^2*a[3]+(4*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]+(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]*a[2]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+(5*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^4*alpha[0]*a[4]+(4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*lambda*a[5]*alpha[0]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]+3*beta[0]^2*alpha[0]*a[2]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2-(1/4)*lambda*beta[0]^2*a[1]-9*mu^2*alpha[1]^2*a[1]*(1/4)+3*mu*a[1]*alpha[0]*beta[0]*(1/2)+(1/4)*(3*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[1]+(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2)*alpha[1]^4*a[3]-w*beta[0]^2-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4] = 0

-k^2*a[1]*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-w*beta[0]^2-(9/4)*mu^2*alpha[1]^2*a[1]+6*beta[0]^2*alpha[0]^2*a[3]-(1/4)*lambda*beta[0]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]+(3/4)*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[1]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+10*beta[0]^2*alpha[0]^3*a[4]+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0] = 0

(1)

numer(lhs(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+(3/4)*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[1]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]+3*beta[0]^2*alpha[0]*a[2]-(9/4)*mu^2*alpha[1]^2*a[1]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2-(1/4)*lambda*beta[0]^2*a[1] = 0))*denom(rhs(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+(3/4)*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[1]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]+3*beta[0]^2*alpha[0]*a[2]-(9/4)*mu^2*alpha[1]^2*a[1]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2-(1/4)*lambda*beta[0]^2*a[1] = 0)) = numer(rhs(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+(3/4)*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[1]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]+3*beta[0]^2*alpha[0]*a[2]-(9/4)*mu^2*alpha[1]^2*a[1]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2-(1/4)*lambda*beta[0]^2*a[1] = 0))*denom(lhs(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+(3/4)*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[1]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]+3*beta[0]^2*alpha[0]*a[2]-(9/4)*mu^2*alpha[1]^2*a[1]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2-(1/4)*lambda*beta[0]^2*a[1] = 0))

-40*lambda^3*B[1]^2*a[4]*alpha[0]*alpha[1]^4+40*lambda^3*B[2]^2*a[4]*alpha[0]*alpha[1]^4-8*lambda^3*B[1]^2*a[3]*alpha[1]^4+8*lambda^3*B[2]^2*a[3]*alpha[1]^4+40*lambda^2*B[1]^2*a[4]*alpha[0]^3*alpha[1]^2-40*lambda^2*B[2]^2*a[4]*alpha[0]^3*alpha[1]^2-4*k^2*lambda^2*B[1]^2*a[1]*alpha[1]^2+4*k^2*lambda^2*B[2]^2*a[1]*alpha[1]^2-16*lambda^3*B[1]^2*a[5]*alpha[0]*alpha[1]^2+16*lambda^3*B[2]^2*a[5]*alpha[0]*alpha[1]^2-80*lambda^2*mu*a[4]*alpha[1]^4*beta[0]+24*lambda^2*B[1]^2*a[3]*alpha[0]^2*alpha[1]^2-24*lambda^2*B[2]^2*a[3]*alpha[0]^2*alpha[1]^2+120*lambda*mu^2*a[4]*alpha[0]*alpha[1]^4-4*lambda^3*B[1]^2*a[1]*alpha[1]^2+4*lambda^3*B[2]^2*a[1]*alpha[1]^2+12*lambda^2*B[1]^2*a[2]*alpha[0]*alpha[1]^2-12*lambda^2*B[2]^2*a[2]*alpha[0]*alpha[1]^2-120*lambda^2*a[4]*alpha[0]*alpha[1]^2*beta[0]^2+24*lambda*mu^2*a[3]*alpha[1]^4+240*lambda*mu*a[4]*alpha[0]^2*alpha[1]^2*beta[0]-40*mu^2*a[4]*alpha[0]^3*alpha[1]^2+4*k^2*mu^2*a[1]*alpha[1]^2-28*lambda^2*mu*a[5]*alpha[1]^2*beta[0]-4*lambda^2*w*B[1]^2*alpha[1]^2+4*lambda^2*w*B[2]^2*alpha[1]^2-24*lambda^2*a[3]*alpha[1]^2*beta[0]^2+32*lambda*mu^2*a[5]*alpha[0]*alpha[1]^2+96*lambda*mu*a[3]*alpha[0]*alpha[1]^2*beta[0]+40*lambda*a[4]*alpha[0]^3*beta[0]^2-24*mu^2*a[3]*alpha[0]^2*alpha[1]^2-4*k^2*lambda*a[1]*beta[0]^2-8*lambda^2*a[5]*alpha[0]*beta[0]^2+7*lambda*mu^2*a[1]*alpha[1]^2+24*lambda*mu*a[2]*alpha[1]^2*beta[0]+12*lambda*mu*a[5]*alpha[0]^2*beta[0]+24*lambda*a[3]*alpha[0]^2*beta[0]^2-12*mu^2*a[2]*alpha[0]*alpha[1]^2-lambda^2*a[1]*beta[0]^2+6*lambda*mu*a[1]*alpha[0]*beta[0]+12*lambda*a[2]*alpha[0]*beta[0]^2+4*mu^2*w*alpha[1]^2-4*lambda*w*beta[0]^2 = 0

(2)

simplify(-40*lambda^3*B[1]^2*a[4]*alpha[0]*alpha[1]^4+40*lambda^3*B[2]^2*a[4]*alpha[0]*alpha[1]^4-8*lambda^3*B[1]^2*a[3]*alpha[1]^4+8*lambda^3*B[2]^2*a[3]*alpha[1]^4+40*lambda^2*B[1]^2*a[4]*alpha[0]^3*alpha[1]^2-40*lambda^2*B[2]^2*a[4]*alpha[0]^3*alpha[1]^2-4*k^2*lambda^2*B[1]^2*a[1]*alpha[1]^2+4*k^2*lambda^2*B[2]^2*a[1]*alpha[1]^2-16*lambda^3*B[1]^2*a[5]*alpha[0]*alpha[1]^2+16*lambda^3*B[2]^2*a[5]*alpha[0]*alpha[1]^2-80*lambda^2*mu*a[4]*alpha[1]^4*beta[0]+24*lambda^2*B[1]^2*a[3]*alpha[0]^2*alpha[1]^2-24*lambda^2*B[2]^2*a[3]*alpha[0]^2*alpha[1]^2+120*lambda*mu^2*a[4]*alpha[0]*alpha[1]^4-4*lambda^3*B[1]^2*a[1]*alpha[1]^2+4*lambda^3*B[2]^2*a[1]*alpha[1]^2+12*lambda^2*B[1]^2*a[2]*alpha[0]*alpha[1]^2-12*lambda^2*B[2]^2*a[2]*alpha[0]*alpha[1]^2-120*lambda^2*a[4]*alpha[0]*alpha[1]^2*beta[0]^2+24*lambda*mu^2*a[3]*alpha[1]^4+240*lambda*mu*a[4]*alpha[0]^2*alpha[1]^2*beta[0]-40*mu^2*a[4]*alpha[0]^3*alpha[1]^2+4*k^2*mu^2*a[1]*alpha[1]^2-28*lambda^2*mu*a[5]*alpha[1]^2*beta[0]-4*lambda^2*w*B[1]^2*alpha[1]^2+4*lambda^2*w*B[2]^2*alpha[1]^2-24*lambda^2*a[3]*alpha[1]^2*beta[0]^2+32*lambda*mu^2*a[5]*alpha[0]*alpha[1]^2+96*lambda*mu*a[3]*alpha[0]*alpha[1]^2*beta[0]+40*lambda*a[4]*alpha[0]^3*beta[0]^2-24*mu^2*a[3]*alpha[0]^2*alpha[1]^2-4*k^2*lambda*a[1]*beta[0]^2-8*lambda^2*a[5]*alpha[0]*beta[0]^2+7*lambda*mu^2*a[1]*alpha[1]^2+24*lambda*mu*a[2]*alpha[1]^2*beta[0]+12*lambda*mu*a[5]*alpha[0]^2*beta[0]+24*lambda*a[3]*alpha[0]^2*beta[0]^2-12*mu^2*a[2]*alpha[0]*alpha[1]^2-lambda^2*a[1]*beta[0]^2+6*lambda*mu*a[1]*alpha[0]*beta[0]+12*lambda*a[2]*alpha[0]*beta[0]^2+4*mu^2*w*alpha[1]^2-4*lambda*w*beta[0]^2 = 0, 'symbolic')

-40*(B[1]-B[2])*((a[4]*alpha[0]+(1/5)*a[3])*alpha[1]^2+(2/5)*a[5]*alpha[0]+(1/10)*a[1])*alpha[1]^2*(B[1]+B[2])*lambda^3+4*(-20*a[4]*beta[0]*alpha[1]^4*mu+(10*(B[1]^2-B[2]^2)*a[4]*alpha[0]^3+6*a[3]*(B[1]^2-B[2]^2)*alpha[0]^2+3*(B[1]^2*a[2]-B[2]^2*a[2]-10*a[4]*beta[0]^2)*alpha[0]-6*beta[0]^2*a[3]-7*a[5]*beta[0]*mu-(B[1]-B[2])*(B[1]+B[2])*(k^2*a[1]+w))*alpha[1]^2-2*(a[5]*alpha[0]+(1/8)*a[1])*beta[0]^2)*lambda^2+(120*(a[4]*alpha[0]+(1/5)*a[3])*mu^2*alpha[1]^4+(240*a[4]*beta[0]*alpha[0]^2*mu+32*(mu^2*a[5]+3*mu*a[3]*beta[0])*alpha[0]+24*beta[0]*mu*a[2]+7*mu^2*a[1])*alpha[1]^2-4*(-10*a[4]*beta[0]*alpha[0]^3+3*(-mu*a[5]-2*a[3]*beta[0])*alpha[0]^2+3*(-beta[0]*a[2]-(1/2)*mu*a[1])*alpha[0]+beta[0]*(k^2*a[1]+w))*beta[0])*lambda+4*alpha[1]^2*mu^2*(-10*a[4]*alpha[0]^3+k^2*a[1]-6*a[3]*alpha[0]^2-3*a[2]*alpha[0]+w) = 0

 

 

 

Error, (in collect) invalid input: collect uses a 2nd argument, x, which is missing

 

Q1 := collect(%, {B__1, B__2})

-40*(B[1]-B[2])*((a[4]*alpha[0]+(1/5)*a[3])*alpha[1]^2+(2/5)*a[5]*alpha[0]+(1/10)*a[1])*alpha[1]^2*(B[1]+B[2])*lambda^3+4*(-20*a[4]*beta[0]*alpha[1]^4*mu+(10*(B[1]^2-B[2]^2)*a[4]*alpha[0]^3+6*a[3]*(B[1]^2-B[2]^2)*alpha[0]^2+3*(B[1]^2*a[2]-B[2]^2*a[2]-10*a[4]*beta[0]^2)*alpha[0]-6*beta[0]^2*a[3]-7*a[5]*beta[0]*mu-(B[1]-B[2])*(B[1]+B[2])*(k^2*a[1]+w))*alpha[1]^2-2*(a[5]*alpha[0]+(1/8)*a[1])*beta[0]^2)*lambda^2+(120*(a[4]*alpha[0]+(1/5)*a[3])*mu^2*alpha[1]^4+(240*a[4]*beta[0]*alpha[0]^2*mu+32*(mu^2*a[5]+3*mu*a[3]*beta[0])*alpha[0]+24*beta[0]*mu*a[2]+7*mu^2*a[1])*alpha[1]^2-4*(-10*a[4]*beta[0]*alpha[0]^3+3*(-mu*a[5]-2*a[3]*beta[0])*alpha[0]^2+3*(-beta[0]*a[2]-(1/2)*mu*a[1])*alpha[0]+beta[0]*(k^2*a[1]+w))*beta[0])*lambda+4*alpha[1]^2*mu^2*(-10*a[4]*alpha[0]^3+k^2*a[1]-6*a[3]*alpha[0]^2-3*a[2]*alpha[0]+w) = 0

(3)

latex(Q1)

-40 \left(B_{1}-B_{2}\right) \left(\left(a_{4} \alpha_{0}+\frac{a_{3}}{5}\right) \alpha_{1}^{2}+\frac{2 a_{5} \alpha_{0}}{5}+\frac{a_{1}}{10}\right) \alpha_{1}^{2} \left(B_{1}+B_{2}\right) \lambda^{3}+4 \left(-20 a_{4} \beta_{0} \alpha_{1}^{4} \mu +\left(10 \left(B_{1}^{2}-B_{2}^{2}\right) a_{4} \alpha_{0}^{3}+6 a_{3} \left(B_{1}^{2}-B_{2}^{2}\right) \alpha_{0}^{2}+3 \left(B_{1}^{2} a_{2}-B_{2}^{2} a_{2}-10 a_{4} \beta_{0}^{2}\right) \alpha_{0}-6 \beta_{0}^{2} a_{3}-7 a_{5} \beta_{0} \mu -\left(B_{1}-B_{2}\right) \left(B_{1}+B_{2}\right) \left(k^{2} a_{1}+w \right)\right) \alpha_{1}^{2}-2 \left(a_{5} \alpha_{0}+\frac{a_{1}}{8}\right) \beta_{0}^{2}\right) \lambda^{2}+\left(120 \left(a_{4} \alpha_{0}+\frac{a_{3}}{5}\right) \mu^{2} \alpha_{1}^{4}+\left(240 a_{4} \beta_{0} \alpha_{0}^{2} \mu +32 \left(\mu^{2} a_{5}+3 \mu  a_{3} \beta_{0}\right) \alpha_{0}+24 \beta_{0} \mu  a_{2}+7 \mu^{2} a_{1}\right) \alpha_{1}^{2}-4 \left(-10 a_{4} \beta_{0} \alpha_{0}^{3}+3 \left(-\mu  a_{5}-2 a_{3} \beta_{0}\right) \alpha_{0}^{2}+3 \left(-\beta_{0} a_{2}-\frac{\mu  a_{1}}{2}\right) \alpha_{0}+\beta_{0} \left(k^{2} a_{1}+w \right)\right) \beta_{0}\right) \lambda +4 \alpha_{1}^{2} \mu^{2} \left(-10 a_{4} \alpha_{0}^{3}+k^{2} a_{1}-6 a_{3} \alpha_{0}^{2}-3 a_{2} \alpha_{0}+w \right)
 = 0

 
 

NULL

Download coment.mw

Major deficiency in Physics[Vectors]; Distinct sets of basis vectors are not recognized!

You can't define vectors in alternative bases like: {\hat{i}',\hat{j}',\hat{k}'} or {\hat{i}_{1},\hat{j}_{2},\hat{k}_{3}}.

This deficiency has been around for a while. I have found other posts regarding this problem.

The deficiency greatly reduces the allowable calculations with Physics[Vector].

Are there any plans to fix this?

Here is my example which shows this deficiency in more detail.

physics_vectors_and_multiple_unit_vectors.mw
 

restart

NULL

NULL

with(Physics[Vectors])

[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface, ParametrizeVolume, Setup, Simplify, `^`, diff, int]

(1)

NULL

Crucial Deficiency in Physics[Vectors]

 

NULL

I can only guess the purpose of the Physics[Vectors] package from reviewing it's corresponding help documentation. My interpretation of the documentation leads me to believe that the package is best used for generating vector equation formulas in different coordinate bases of a SINGLE coordinate system.

 

This means one can easily generate position vector expressions such as:

 

r_ = _i*x+_j*y+_k*z

r_ = _i*x+_j*y+_k*z

(1.1)

Cylindrical Position Vector

 

The position vector in a cylindrical basis is given by:

 

r_ = ChangeBasis(rhs(r_ = _i*x+_j*y+_k*z), 2)

r_ = (x*cos(phi)+y*sin(phi))*_rho+(cos(phi)*y-sin(phi)*x)*_phi+z*_k

(1.1.1)

r_ = ChangeBasis(rhs(r_ = _i*x+_j*y+_k*z), 2, alsocomponents)

r_ = _k*z+_rho*rho

(1.1.2)

NULL

NULLNULLNULL

Spherical Position Vector

 

NULL

r_ = ChangeBasis(rhs(r_ = _i*x+_j*y+_k*z), 3)

r_ = (y*sin(phi)*sin(theta)+x*sin(theta)*cos(phi)+z*cos(theta))*_r+(y*sin(phi)*cos(theta)+x*cos(phi)*cos(theta)-z*sin(theta))*_theta+(cos(phi)*y-sin(phi)*x)*_phi

(1.2.1)

r_ = ChangeBasis(rhs(r_ = _i*x+_j*y+_k*z), 3, alsocomponents)

r_ = r*_r

(1.2.2)

NULL

NULL

As is known from the vector analysis of curvilinear coordinate systems the basis vectors can depend on the coordinates in question.

 

In cylindrical, the basis vectors are

 

_rho = ChangeBasis(_rho, 1)

_rho = _i*cos(phi)+sin(phi)*_j

(1.2)

_phi = ChangeBasis(_phi, 1)

_phi = -sin(phi)*_i+cos(phi)*_j

(1.3)

and in spherical, the basis vectors are

 

_r = ChangeBasis(_r, 1)

_r = sin(theta)*cos(phi)*_i+sin(theta)*sin(phi)*_j+cos(theta)*_k

(1.4)

_theta = ChangeBasis(_theta, 1)

_theta = cos(theta)*cos(phi)*_i+cos(theta)*sin(phi)*_j-sin(theta)*_k

(1.5)

_phi = ChangeBasis(_phi, 1)

_phi = -sin(phi)*_i+cos(phi)*_j

(1.6)

NULL

NULL

NULL

Example of this Deficiency using Biot-Savart Law

 

NULL

Biot-Savart law can be used to calculate a magnetic field due to a current carrying wire. The deficiency in question can be observed by explicity constructing the integrand in the Biot-Savart integral defined below.

NULL

NULL

NULL

In electrodynamics, quantum mechanics and applied mathematics, it is common practice to define a position of observation by a vector `#mover(mi("r"),mo("&rarr;"))` and a position of the source responsible for generating the field by a vector diff(`#mover(mi("r"),mo("&rarr;"))`(x), x).

 

It is just as common to define the difference in these vectors as

 

l_ = r_-(diff(r(x), x))*_

l_ = r_-`r'_`

(1.3.1)

and thus

 

dl_ = dr_-(diff(dr(x), x))*_

dl_ = dr_-`dr'_`

(1.3.2)

as found in the integrand of the Biot-Savart integral.

NULL

It suffices to consider `#mover(mi("l"),mo("&rarr;"))` = `#mover(mi("r"),mo("&rarr;"))`-`#mover(mi("r'"),mo("&rarr;"))` in a cylindrical basis for this argument.

 

The observation position is:

 

`#mover(mi("r"),mo("&rarr;"))` = rho*`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`+z*`#mover(mi("k"),mo("&and;"))`

NULL

The source position is:

 

diff(`#mover(mi("r"),mo("&rarr;"))`(x), x) = (diff(z(x), x))*(diff(`#mover(mi("k"),mo("&and;"))`(x), x))+(diff(rho(x), x))*(diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x))

NULL

`#mover(mi("l"),mo("&rarr;"))` = `#mover(mi("r"),mo("&rarr;"))`-(diff(`#mover(mi("r"),mo("&rarr;"))`(x), x)) and `#mover(mi("r"),mo("&rarr;"))`-(diff(`#mover(mi("r"),mo("&rarr;"))`(x), x)) = z(x)*`#mover(mi("k"),mo("&and;"))`-(diff(z(x), x))*(diff(`#mover(mi("k"),mo("&and;"))`(x), x))+rho*`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`-(diff(rho(x), x))*(diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x))

NULL

The deficiency in question arises because MAPLE cannot define multiple unit vectors in distinct bases such as {`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`, diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x)} or {`#mscripts(mi("&rho;",fontstyle = "normal"),mn("1"),none(),none(),mo("&and;"),none(),none())`, `#mscripts(mi("&rho;",fontstyle = "normal"),mn("2"),none(),none(),mo("&and;"),none(),none())`}.  These pairs of unit vectors arise naturally, as shown above in Biot-Savart law.

NULL

If we look at `#mover(mi("&rho;",fontstyle = "normal"),mo("&circ;"))` and  diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&circ;"))`(x), x) generally, they look like:

NULL

`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))` = `#mover(mi("i"),mo("&and;"))`*cos(phi)+sin(phi)*`#mover(mi("j"),mo("&and;"))`

NULL

diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x) = (diff(`#mover(mi("i"),mo("&and;"))`(x), x))*cos(diff(phi(x), x))+sin(diff(phi(x), x))*(diff(`#mover(mi("j"),mo("&and;"))`(x), x))

NULL

If the bases vectors {`#mover(mi("i"),mo("&and;"))`, `#mover(mi("j"),mo("&and;"))`, `#mover(mi("k"),mo("&and;"))`} and {diff(`#mover(mi("i"),mo("&and;"))`(x), x), diff(`#mover(mi("j"),mo("&and;"))`(x), x), diff(`#mover(mi("k"),mo("&and;"))`(x), x)} are Cartesian and are not related related through rotations so that

NULL

"(i)*i' =(|i|)*|i'|*cos(0)=1"``NULL

NULL

"(j)*(j)' =(|j|)*|(j)'|*cos(0)=1"NULL

NULL

"(k)*(k)' =(|k|)*|(k)'|*cos(0)=1 "

NULL

and so,NULL

 

`#mover(mi("i"),mo("&circ;"))` = diff(`#mover(mi("i"),mo("&circ;"))`(x), x)

NULL

`#mover(mi("j"),mo("&circ;"))` = diff(`#mover(mi("j"),mo("&circ;"))`(x), x)

NULL

`#mover(mi("k"),mo("&circ;"))` = diff(`#mover(mi("k"),mo("&circ;"))`(x), x)

NULL

the radial unit vectors in cylindrical are then,

 

`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))` = `#mover(mi("i"),mo("&and;"))`*cos(phi)+sin(phi)*`#mover(mi("j"),mo("&and;"))`

NULL

diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x) = `#mover(mi("i"),mo("&and;"))`*cos(diff(phi(x), x))+sin(diff(phi(x), x))*`#mover(mi("j"),mo("&and;"))`

NULL

In typical problems, the anglular location of the observation point, φ, is distinct from the angular location of the source, diff(phi(x), x) and so under this condition, `#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))` <> diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x).

 

Consider the classic problem of the magnetic field due to a circular current carrying wire. Surely, the angular coordinate of one location of the current carrying wire  is different from the angular coordinate  of an observation point hovering above and off-axis on the other side of the current carrying wire. See figure below.

NULL

NULL

NULL

NULL

Therefore,

 

`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))` <> diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x)

NULL

NULL

What happens in MAPLE when you try to define a second distinct unit vector diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x)?

NULL

One can easily find `#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`.

NULL

_rho

_rho

(1.3.3)

NULL

NULLIf you try to define diff(`#mover(mi("&rho;",fontstyle = "normal"),mo("&and;"))`(x), x) ...

 

 

diff(_rho(x), x)

`_rho'`

(1.3.4)

So using a prime doesn't work.

NULL

You could try a numbered subscript...

`_&rho;__2`

_rho__2

(1.3.5)

but that doesn't work.

 

You could try an indexed unit vector...

NULL

_rho[2]

_rho[2]

(1.3.6)

which can be define but is not recognized by Physics[Vectors] since...

 

NULL

ChangeBasis(_rho[2], 1)

Error, (in Physics:-Vectors:-Identify) incorrect indexed use of a unit vector: _rho[2]

 

NULL

And so it's just not possible with the current implementation.

``

``

NULL

NULL


 

Download physics_vectors_and_multiple_unit_vectors.mw

 

 

Where can I found details about Statistics:-Sample(..., method=envelope).

It would be nice to have a link to a description of the envelope method Sample uses.
For instance does it share some features of the Cuba library for numeric integration? Does it use the same envelope method evalf/Int(..., method=_CubaSuave)) uses?

Thanks in advance.

This is my first time solving an equation where the solution includes a function W(t). I'm not sure how to work with it. Does anyone have insights about this function? I also don't have much knowledge about stochastic processes, which I think might be related. How can I gather enough information to understand and plot such a function?
in that paper i saw he talk about wiener process!

restart;

local gamma;

gamma

(1)

``

T3 := (B[1]*(coth(2*n^2*(delta^2-w)*k*t/((k*n-1)*(k*n+1))+x)-1))^(1/(2*n))*exp(I*(-k*x+w*t+delta*W(t)-delta^2*t))

(B[1]*(coth(2*n^2*(delta^2-w)*k*t/((k*n-1)*(k*n+1))+x)-1))^((1/2)/n)*exp(I*(-k*x+w*t+delta*W(t)-delta^2*t))

(2)

``

params := {B[1]=1,n=2,delta=1,w=1,k=3 };

{delta = 1, k = 3, n = 2, w = 1, B[1] = 1}

(3)

``

insert numerical values

solnum :=subs(params, T3);

(coth(x)-1)^(1/4)*exp(I*(-3*x+W(t)))

(4)

``

P := Array(1 .. 2); P[1] := plot3d(map(Re, solnum), x = -2 .. 2, t = 0 .. 10, title = Re); P[2] := plot3d(map(Im, solnum), x = -20 .. 20, t = -10 .. 10, title = Im); plots:-display(P)

 

 

 

 

 

plot3d(map(abs, solnum), x = -10 .. 10, t = 0 .. 5)

 

 

Download graph-stochastic.mw

I'm having trouble solving this system of differential equations. I haven't solved systems of differential equations before but i tried defining the system and then using dsolve, but it couldn't solve all the equations.

Hope you can help.

NULL

diff(Q1(t), t) = -k1*Q1(t)

 

diff(Q2(t), t) = k1*Q1+k3/Q2(t)-k2*Q2(t)-k4*Q2(t)

 

diff(Q3(t), t) = k4*Q2

 

diff(Q4(t), t) = k2*Q2-k3/Q2

 

NULL

Download System_Of_Differential_Equations.mw

I'm trying to transform a partial differential equation (PDE) into an ordinary differential equation (ODE) as demonstrated in the paper. However, I find some steps confusing and difficult to follow. The process often feels chaotic, and managing the complexity of the equations is overwhelming. Could you suggest an effective and systematic method to handle such transformations more easily?

restart

with(PDEtools)

with(LinearAlgebra)

NULL

with(SolveTools)

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

(1)

declare(Omega(x, t)); declare(U(xi))

Omega(x, t)*`will now be displayed as`*Omega

 

U(xi)*`will now be displayed as`*U

(2)

tr := {t = tau, x = tau*c[0]+xi, Omega(x, t) = U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))}

{t = tau, x = tau*c[0]+xi, Omega(x, t) = U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))}

(3)

P1 := diff(Omega(x, t)^m, t)

Omega(x, t)^m*m*(diff(Omega(x, t), t))/Omega(x, t)

(4)

L1 := PDEtools:-dchange(tr, P1, [xi, tau, U])

(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*(-((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))*c[0]+I*U(xi)*(-k*c[0]+w+delta*(diff(W(tau), tau))-delta^2)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))/(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))

(5)
 

pde1 := I*(diff(Omega(x, t)^m, t))+alpha*(diff(Omega(x, t)^m, `$`(x, 2)))+I*beta*(diff(abs(Omega(x, t))^(2*n)*Omega(x, t)^m, x))+m*sigma*Omega(x, t)^m*(diff(W(t), t)) = I*gamma*abs(Omega(x, t))^(2*n)*(diff(Omega(x, t)^m, x))+delta*abs(Omega(x, t))^(4*n)*Omega(x, t)^m

I*Omega(x, t)^m*m*(diff(Omega(x, t), t))/Omega(x, t)+alpha*(Omega(x, t)^m*m^2*(diff(Omega(x, t), x))^2/Omega(x, t)^2+Omega(x, t)^m*m*(diff(diff(Omega(x, t), x), x))/Omega(x, t)-Omega(x, t)^m*m*(diff(Omega(x, t), x))^2/Omega(x, t)^2)+I*beta*(2*abs(Omega(x, t))^(2*n)*n*(diff(Omega(x, t), x))*abs(1, Omega(x, t))*Omega(x, t)^m/abs(Omega(x, t))+abs(Omega(x, t))^(2*n)*Omega(x, t)^m*m*(diff(Omega(x, t), x))/Omega(x, t))+m*sigma*Omega(x, t)^m*(diff(W(t), t)) = I*gamma*abs(Omega(x, t))^(2*n)*Omega(x, t)^m*m*(diff(Omega(x, t), x))/Omega(x, t)+delta*abs(Omega(x, t))^(4*n)*Omega(x, t)^m

(6)

NULL

L1 := PDEtools:-dchange(tr, pde1, [xi, tau, U])

I*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*(-((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))*c[0]+I*U(xi)*(-k*c[0]+w+delta*(diff(W(tau), tau))-delta^2)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))/(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))+alpha*((U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m^2*((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^2/(U(xi)^2*(exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^2)+(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*((diff(diff(U(xi), xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-(2*I)*(diff(U(xi), xi))*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-U(xi)*k^2*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))/(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))-(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^2/(U(xi)^2*(exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^2))+I*beta*(2*(abs(U(xi))*exp(-Im(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^(2*n)*n*((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))*abs(1, U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m/(abs(U(xi))*exp(-Im(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))+(abs(U(xi))*exp(-Im(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^(2*n)*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))/(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))))+m*sigma*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*(diff(W(tau), tau)) = I*gamma*(abs(U(xi))*exp(-Im(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^(2*n)*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m*m*((diff(U(xi), xi))*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau))-I*U(xi)*k*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))/(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))+delta*(abs(U(xi))*exp(-Im(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^(4*n)*(U(xi)*exp(I*(-k*(tau*c[0]+xi)+w*tau+delta*W(tau)-delta^2*tau)))^m

(7)

``

``

(8)

Download transform-pde-to-ode-hard_example.mw

I have expected the opposite. Is exp already optimised that hardware floats do not make sense or does the conversion of the argument to hardware floast eats up all the benefit of using hardware floats?

restart;
CodeTools:-Usage( for i from 1 to 100 by 0.1 do exp(i) end do):
CodeTools:-Usage( for i from 1 to 100 by 0.1 do (evalhf@exp)(i) end do):
memory used=1.54MiB, alloc change=0 bytes, cpu time=16.00ms, real time=15.00ms, gc time=0ns
memory used=5.88MiB, alloc change=32.00MiB, cpu time=31.00ms, real time=33.00ms, gc time=0ns
1 2 3 4 5 6 7 Last Page 1 of 2167