Carl Love

Carl Love

28160 Reputation

25 Badges

14 years, 8 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

This was trickier than I expected, probably because of a bug in Maple.

a,b,c:= 'RootOf(x^3-x^2-x-1, x, index= k)' $ k= 1..3:
evala(
     evala((a^2014-b^2014)/(a-b)) +
     evala((b^2014-c^2014)/(b-c)) +
     evala((c^2014-a^2014)/(c-a))
);

8334826359539004144901989461804527198477505526542440583547585702\
  52333926756888184793278322997734194204312479082099178141286346\
  19093619632692976466847552803713809622682218736343227666160495\
  08221814701281712551536720818416541447585602970496630535429027\
  14035268777406803186312918798664076365740099763393011180604524\
  46886004153088055743393327297340352231807305287646640702341995\
  05000473477741234080182885893584028536441114628416893583917671\
  56862792958972093429609456884388481631837778020024880944176069\
  54812869711446457823816195550797524

The answer is just a large integer. It is a[2014] for the Fibonacci-like recurrence a[n] = a[n-1]+a[n-2]+a[n-3], a[0]=0, a[1]=3, a[2]=2.

The bug ensues if you try to do it with a single evala command. It seems to get stuck in an infinite loop.


(**)

restart:

Note that Maple's series command will return the linear and quadratic approximations. (The taylor command does the same thing.)

Linear is considered order 2:

(**)

convert(series(f(x), x=a, 2), polynom);

f(a)+(D(f))(a)*(x-a)

Quadratic is order 3:

(**)

convert(series(f(x), x=a, 3), polynom);

f(a)+(D(f))(a)*(x-a)+(1/2)*((D@@2)(f))(a)*(x-a)^2

Applying this to your function, we get:

(**)

P1:= convert(series(arccos(x), x=0, 2), polynom);

(1/2)*Pi-x

(**)

P2:= convert(series(arccos(x), x=0, 3), polynom);

(1/2)*Pi-x

They are the same because x=0 is an inflection point, which will be clear from the plot.

(**)

plot([arccos(x), P1, P2], x= -1..1);

(**)

 


Download linear_approx.mw


(**)

restart:

(**)

Digits:= 2^6:

(**)

eq:= 6.750969101*10^26*_Z*1000^(9/59)*800^(41/59)-1.015701668*10^27*1000^(50/59)-1.189275761*10^27*800^(50/59)+1.103397224*10^26*_Z*1000^(32/59)*800^(18/59)-7.457053114*10^23*3^(1/5)*1000^(41/59)*800^(50/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-7.457053114*10^23*3^(1/5)*1000^(50/59)*800^(41/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-1.988547497*10^26*3^(1/5)*1000^(32/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-2.485684371*10^26*3^(1/5)*800^(32/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)+5.575397080*10^26*_Z*1000^(41/59)*800^(9/59)+1.308382216*10^26*_Z*1000^(18/59)*800^(32/59)+1.378754484*10^27*_Z*800^(50/59)+1.114796623*10^27*_Z*1000^(50/59)-1.144263788*10^26*1000^(32/59)*800^(18/59)-1.028573111*10^26*1000^(18/59)*800^(32/59)-5.399913635*10^26*1000^(41/59)*800^(9/59)-5.544622178*10^26*1000^(9/59)*800^(41/59):

(**)

indets(eq, {name, name^numeric});

{_Z, _Z^(4/5)}

Use a change of variables to convert to a polynomial:

(**)

eq_simp:= simplify(eval(eq, _Z= z^5), symbolic);

1252302292256880108322573586627.994134196842318012573369627002321*z^5-1113835100487674276481580185899.975871697792466972768109672627414-33206779716982073462031352179419335.13149299607560858653432559958*z^4

(**)

r:= fsolve(eq_simp);

26516.58463160465918042279703097330382230532550324190250597416872

Residual in the simplified equation:

(**)

eval(eq_simp, z= r);

0.

Resisual in the original equation:

(**)

evalf(eval(eq, _Z= r^5));

-0.324e-10

(**)

sol1:= r^5;

13109554349246137246800.44849539580049816642933649202738691745839

For comparison, apply RootFinding:-Analytic to the original equation, using the same bounding box as did Markiyan. (This takes much longer than fsolve.)

(**)

sol2:= RootFinding:-Analytic(eq, _Z= 0-10*I..10^30+100*I);

13109554349246137246800.44849539580049816642933649202738691745852

(**)

sol1-sol2;

-0.13e-39

A plot confirms this root.

(**)

plot(eq, _Z= 0..2*10^22);

(**)

 


Download sensitive_equatio.mw

First, it helps tremendously if you indent your code, like this:

restart:
P := array([[1, 4], [2, 3], [3, 2], [4, 1], [6, 5], [6, 1], [6, 2], [5, 3]]);

j := 1; k := j+2; Fold := 0; Fnew := 1; counter := 0;

for i to 6 do
     j := i; k := j+2;
     while Fnew > Fold do
          Fold := Fnew;
          L[j, k] := evalf((P[k, 1]-P[j, 1])^2+(P[k, 2]-P[j, 2])^2)^(1/2);
          m[j, k] := (P[k, 1]-P[j, 1])/(P[k, 2]-P[j, 2]);
          a := m[j, k];
          b := -1;
          c := P[j, 2]-m[j, k]*P[j, 1];
          for i from j+1 to k-1 do
               e[i] := evalf(abs(a*P[i, 1]+b*P[i, 2]+c)/sqrt(a^2+b^2));
               E[j, k] := sum(e[m], m = j+1 .. k-1)
          end do;
          Fnew := L[j, k]-E[j, k];
          counter := counter +1;
          k := k+1
     end do;
     right[j] := k-2
 end do;

And when I look at it like that, I immediately see the problem: You used the same index variable, i, for both the inner and outer for loops. So I simply changed the loop index to ii for the inner loop, and I changed its three invocations within the loop to ii also:

          for ii from j+1 to k-1 do 
               e[ii] := evalf(abs(a*P[ii, 1]+b*P[ii, 2]+c)/sqrt(a^2+b^2));
               E[j, k] := sum(e[m], m = j+1 .. k-1)
          end do;

And then I got 6 values in right.

I realize that this is your first time. In the future please upload a worksheet instead of using a screen shot. The people who answer questions here generally do not take kindly to having to retype your formulas from a screen shot.

There might be a difference between Maple 13 and 17 at issue here. I did not get the warning, and I got values for Z. It's just a cubic equation. Here are three techniques to solve for Z. Almost surely, at least one of these will work for you.

restart:
Digits:= 10:
P:= [beta= 0.08711686224, q= 4.446850237, epsilon= 1-sqrt(2), sigma= 1+sqrt(2)]:
EC1:= Z=1+beta-(q*beta*(Z-beta)/(Z+epsilon*beta)/(Z+sigma*beta));

Three solution techniques:
1. Solve symbolically first, then substitute parameters:
evalf(eval([solve(EC1, Z)], P));
         [0.6906702280, 0.1111064548 + 0.1567590067 I,
           0.1111064548 - 0.1567590067 I]


2. Substitute parameters first, then solve. This is what you were trying.
solve(eval(EC1, P), Z);
         0.6906702273, 0.1111064551 + 0.1567590075 I,
           0.1111064551 - 0.1567590075 I


3. Use numeric solver fsolve:
fsolve(eval(EC1, P), Z);
                          0.6906702273


The package name is Orbitals, not "orbital". Does the folder C:\MapleLibs contain the files Orbitals.lib, Orbitals.hdb, and Oribitals.ind?

To send the output to a file, include the filename option:

codegen[fortran](A, filename= "myfortranfile.txt");

Try to convert your code to use the newer CodeGeneration:-Fortran, as the older one is deprecated. In that case, to get file output use

CodeGeneration:-Fortran(A, output= "myfortranfile.txt");

I told you in my Answer to your previous Question that functions should be defined via f:= (x,y)-> ...and not via f(x,y):= .... Maple will often let you "get away with" using the latter form, but this time it didn't. I can't explain why it didn't work this time, nor can I explain why Maple seems to totally ignore the command (that seems totally strange). But I do know that if you change it to

pathZ:= (r,theta)-> ...

then it will work.

Just change your command from x=0 to h=0:

taylor(y(x+h), h=0);

To have matrix and vector computations done in single precision, create the Matrices and Vectors with the option datatype= float[4].

Example:

V:= Vector(1000, n-> 1/n^3, datatype= float[4]);
V.V;

                     1.01734306201200031

Compare with the result obtained using datatype= float[8], which is double precision.

                     1.01734306198444679

Note that the results begin to differ at the 9th decimal place.

You're right. I don't think that there's any way to get it to accept infinity. Here's a workaround: Use gamma instead of infinity. When the tutor is done, give the command:

limit(subs(gamma= x, %), x= infinity);

For example,

count:= 0:
for i from 1 to 9 do
     for j from 1 to i do
          count:= count+1
     end do
end do;
count;
                                              45

Of course, this example is trivial, and there are much more efficient ways to perform this particular computation.

interface(showassumed= 0):
assume(nu::real, m>0, k>0, T>0):
Int(nu*4*Pi*(m/2/Pi/k/T)^(3/2)*nu^2*exp(-m*nu^2/2/k/T), nu= 0..infinity);

value(%);

Whenever a symbol is redefined by a package, the original value of the symbol can be accessed by prefixing the symbol with :-. So in this case you'd use :-``.

with(Physics):
``(3)*``(2);
Error, (in GetDifferentiationVariables) differentiation variables for dAlembertian, d_ and D_ are not defined; either use Setup(differentiationvariables =  ...) to set them, or indicate the differentiation variables as a second argument to dAlembertian

:-``(3)*:-``(2);
                            (3) (2)

There is a policy stated in the 6th paragraph at ?name : 

Names that start with an underscore (with the exception of _Env) are used as global variable names by Maple and are effectively reserved for use by library code. Do not use a symbol beginning with an underscore. The empty name `` must never be assigned a value.

It seems that Physics violates this.

 

For the constraints to DirectSearch, you could generate all pairs of the variables like this:

map(`<>`@op, combinat:-choose({a,b,c,d}, 2));

First 347 348 349 350 351 352 353 Last Page 349 of 396