acer

32303 Reputation

29 Badges

19 years, 310 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Yes, different forms of an expression can have different amounts of numeric error (eg, loss of precision).

It is quite common in numeric models of evaluation in scientific computing, and is not specific to Maple.

One way to deal with it is to raise working precision (either using the Digits environment variable, or more locally via the manner of calling the evalf command).

restart;

values := {m = 1400, d = 1/2 * 1.225*0.74, xpr = 27+7/9,
           A = -94.25, B = 4096, C = -2699, D = 131947,
           u0 = 0.1, P_max = 9705, u = 0.5}:

Pa_R := ((((A*D - B*C)*(sqrt(4*A^2*d*u*xpr^2 + 4*A*d*(A*D - B*C)*xpr
                        + (A*D - B*C)^2) + (2*d*xpr + D)*A - C*B))/(2*A^2*d)
          + u*xpr^2)*(-A*D + C*B - sqrt(4*A^2*d*u*xpr^2
                      + 4*A*d*(A*D - B*C)*xpr + (A*D - B*C)^2)))/(2*A*d);

(1/2)*((1/2)*(A*D-B*C)*((4*A^2*d*u*xpr^2+4*A*d*(A*D-B*C)*xpr+(A*D-B*C)^2)^(1/2)+(2*d*xpr+D)*A-B*C)/(A^2*d)+u*xpr^2)*(-A*D+B*C-(4*A^2*d*u*xpr^2+4*A*d*(A*D-B*C)*xpr+(A*D-B*C)^2)^(1/2))/(A*d)

Pa_R2 := simplify(Pa_R);

-(1/4)*(A*D-B*C+(4*A^2*d*u*xpr^2+4*A*d*(A*D-B*C)*xpr+(A*D-B*C)^2)^(1/2))*((A*D-B*C)*(4*A^2*d*u*xpr^2+4*A*d*(A*D-B*C)*xpr+(A*D-B*C)^2)^(1/2)+(2*d*u*xpr^2+2*D*d*xpr+D^2)*A^2-2*B*C*(d*xpr+D)*A+B^2*C^2)/(A^3*d^2)

eval(Pa_R, values);

9717.251005

eval(Pa_R2, values);

9314.704850

evalf[15](eval(Pa_R, values));

9717.30008417060

evalf[15](eval(Pa_R2, values));

9717.30163215298

for dd from 10 to 15 do
  forget(evalf):
  printf("   %a, %a\n", dd, evalf[dd](eval(Pa_R, values)));
end do;

   10, 9717.251005
   11, 9722.0166015
   12, 9717.25243950
   13, 9717.347726855
   14, 9717.3048483735
   15, 9717.30008417060

for dd from 10 to 15 do
  forget(evalf):
  printf("   %a, %a\n", dd, evalf[dd](eval(Pa_R2, values)));
end do;

   10, 9314.704850
   11, 9763.1922855
   12, 9714.89360302
   13, 9717.308535045
   14, 9717.3085321810
   15, 9717.30163215298

 

Download num_examp.mw

There are ways to gauge the amount of accrued numeric error. For example the precision of those supplied values for the parameters could be quantified (as tolerances, or intervals, etc), which could allow one to gauge accuracy of the resulting evaluation (as a function of the variable working precision).

Are you trying to multiply both numerator and denominator by the conjugate of the denominator?

Perhaps you also want to follow that up with simplification, and to illustrate doing it step by step?

Is the name x supposed to represent an unknown purely real quantity?

restart;

interface(typesetting=extended):

 

expr := (8*x+3*I)/(5*x-I);

(8*x+3*I)/(5*x-I)

n := numer(expr)/sign(numer(expr));

8*x+3*I

d := denom(expr)/sign(numer(expr));

5*x-I

f := conjugate(d) assuming x::real;

5*x+I

(n %* f) / (d %* f);

`%*`(8*x+3*I, 5*x+I)/`%*`(5*x-I, 5*x+I)

expand(n*f) / expand(d*f);

(40*x^2+(23*I)*x-3)/(25*x^2+1)


Alternatively,

rationalize(expr);

(8*x+3*I)*(5*x+I)/(25*x^2+1)

simplify(expand(%));

(40*x^2+(23*I)*x-3)/(25*x^2+1)

simplify(evalc(expr));

(40*x^2+(23*I)*x-3)/(25*x^2+1)

Download rtnz.mw

See the CodeTools:-Usage or the time commands.

With the time command, you may wish also to utilize assign, eg,

    time[real]( assign('ans', somecomputation) );

It depends on whether you want the timing to be merely printed, or accessible programmatically as a returned value.

You could assign the string (full and explicit name of the .mla file, with path) to a Maple name. Then you could pass that very same name to both commands, and in doing so avoid mismatch mistakes.

For example,

filenm := "C:\\Home\\Maple projects\\Test_package\\test.mla":
LibraryTools:-Create(filenm);
LibraryTools:-Save('testSavelib', filenm);

The above correction would also work in your original code, and avoid your mismatch between the your two lines related to the extra whitespace (that the second line had, while the first did not).

Regardless of the nature of your particular problem, I think that using LibraryTools:-Save is generally better (less error prone) than using savelib.

Why aren't you using fsolve for this?

CodeTools:-Usage(fsolve(1 + 0.15 = (1 + r__x)^1.28858));
memory used=1.33MiB, alloc change=0 bytes, cpu time=16.00ms, real time=16.00ms, gc time=0ns

              0.1145625357

Is there a numeric value for parameter a?

From your grainy image it looks like you may have accidentally only made the equation,

   a = 1

instead of doing the assignment,

   a := 1

It's more useful to upload and attach your actual worksheet (green up-arrow in the Mapleprimes editor) than it is to show only an image of your code.

Are you hoping to utilize all of those equations rC1,rC2,rf1,rA,rh to substitute for all of C1,C2,f1,A, and h?

restart

with(plots):

``

psi := I*(C1*exp(A)-C2*exp(-A))*exp(-I*t*(1/2))/f1;

I*(C1*exp(A)-C2*exp(-A))*exp(-((1/2)*I)*t)/f1

rf1 := f1 = sqrt(h^2-1);

f1 = (h^2-1)^(1/2)

rh := h = f^2+1;

h = f^2+1

rA := `assuming`([A = simplify(subs(rh, sqrt(h^2-1)*(x+I*h*t)))], [f > 0]);

A = f*(f^2+2)^(1/2)*(I*f^2*t+I*t+x)

rC1 := `assuming`([C1 = simplify(subs(rf1, rh, sqrt(h-f1)))], [f > 0]);

C1 = (f^2+1-f*(f^2+2)^(1/2))^(1/2)

rC2 := `assuming`([C2 = simplify(subs(rf1, rh, sqrt(h+f1)))], [f > 0]);

C2 = (f^2+1+f*(f^2+2)^(1/2))^(1/2)

psi_f := eval[recurse](psi, [rC1, rC2, rf1, rA, rh]);

I*((f^2+1-f*(f^2+2)^(1/2))^(1/2)*exp(f*(f^2+2)^(1/2)*(I*f^2*t+I*t+x))-(f^2+1+f*(f^2+2)^(1/2))^(1/2)*exp(-f*(f^2+2)^(1/2)*(I*f^2*t+I*t+x)))*exp(-((1/2)*I)*t)/((f^2+1)^2-1)^(1/2)

simplify(limit(psi_f, f = 0, right));

exp(-((1/2)*I)*t)*((2*I)*x-I-2*t)

simplify(limit(psi_f, f = 0, left));

-exp(-((1/2)*I)*t)*((2*I)*x-I-2*t)

 

NULL

Download limit1.mw

At lower (and default) working precision this example can encounter enough numeric round-off error that the symbolic solve command finds spurious roots if you perform non-inert symbolic solving and integration on the expression containing floating-point coefficients.

I suggest:
1) don't use RealDomain (it's not doing what you think, anyway, in your example).
2) use fsolve for numeric root-finding instead of solve here.
3) Don't introduce floating-point values instead of exact rationals if you can practically avoid doing so, especially if you intend on attempting symbolic integration or solving. (Symbolic integration/solving of an expression containing floats can sometimes do worse than trying to solve it numerically, or trying to solve an exact problem.)

If you left out RealDomain, and used exact rationals, then the spurious solutions from solve (due to floating-point round-off error) could be avoided by raising Digits. But it'd be better and faster to simply use a numeric root-finder for this problem which likely doesn't even have an explicit symbolic solution.

restart

f := proc (x) options operator, arrow; 6-sqrt(-x^2+8*x+9) end proc

eq := 10 = int(f(x), x = 0 .. L)

10 = -6-(25/2)*arcsin(4/5)+6*L-(1/2)*(-L^2+8*L+9)^(1/2)*L+2*(-L^2+8*L+9)^(1/2)-(25/2)*arcsin((1/5)*L-4/5)

fsolve(eq, L)

6.810993081

fsolve(eq, L = -10 .. 10, maxsols = 5)

6.810993081

Download Working_ac.mw

 

If you really want to see more of the specifics of how it went wrong...

restart

f := proc (x) options operator, arrow; 6-sqrt(-x^2+8*x+9) end proc

proc (x) options operator, arrow; 6-sqrt(-x^2+8*x+9) end proc

eq := 10 = int(f(x), x = 0 .. L)

10 = -6-(25/2)*arcsin(4/5)+6*L-(1/2)*(-L^2+8*L+9)^(1/2)*L+2*(-L^2+8*L+9)^(1/2)-(25/2)*arcsin((1/5)*L-4/5)

oof := 10 = int(evalf(f(x)), x = 0 .. L)

10 = -17.59119023+6.*L-.5000000000*(-1.*L^2+8.*L+9.)^(1/2)*L+2.*(-1.*L^2+8.*L+9.)^(1/2)-12.50000000*arcsin(.2000000000*L-.8000000000)

feq := evalf(eq)

10. = -17.59119022+6.*L-.5000000000*(-1.*L^2+8.*L+9.)^(1/2)*L+2.*(-1.*L^2+8.*L+9.)^(1/2)-12.50000000*arcsin(.2000000000*L-.8000000000)

solve(eq, L)

5*sin(RootOf(25*(cos(_Z)^2)^(1/2)*sin(_Z)-60*sin(_Z)+25*_Z+25*arcsin(4/5)-16))+4

RealDomain:-solve(feq, L)

8.855609656, -.4725737212, 6.810993078

solve(feq, L)

8.855609656, -.4725737212, 6.810993078, -.7764912837+15.69863487*I, 8.138450513-15.36358890*I, -.7764912837-15.69863487*I, 8.138450513+15.36358890*I, 2.341723694-6.481846912*I, 2.341723694+6.481846912*I, 13.34312523-18.81939387*I, -5.700438002+19.09985409*I, 13.34312523+18.81939387*I, -5.700438002-19.09985409*I, 16.71699261-17.45692870*I, -9.003420055+17.75423373*I, 16.71699261+17.45692870*I, -9.003420055-17.75423373*I

If we approximate eq at higher working precision then terms like arcsin(4/5) don't incur as must loss of accuracy (which propagates forward as greater inaccuracy, leading to spurious solutions).
Digits := 30; solve(evalf(eq))

6.81099308113251113676562942220

NULL

Download Not_working_ac.mw

You can call fsolve to compute a root between each of those.

But don't try and do it by appending (to a list, or Vector). That's a poor way to program, and habitually leads to unnecessary inefficiency in programs as they scale.

Instead, you could simply use seq (or map) to form the list. Or you could use Vector with a constructor operator.

restart

gamma1 := .1093:

f := q-(1-alpha)*tan(q)-c*tan(q)

q-0.4646295e-3*tan(q)-0.2298931352e-4*tan(q)/(0.3800000000e-3*q^2-0.2300000000e-4)

n := 10:

Asymptotes

fOdd := proc (j) options operator, arrow; (1/2)*(2*j+1)*Pi end proc; OddAsymptotes := Vector[row](n, fOdd)

fOdd := proc (j) options operator, arrow; (j+1/2)*Pi end proc

OddAsymptotes := Vector[row](10, {(1) = (3/2)*Pi, (2) = (5/2)*Pi, (3) = (7/2)*Pi, (4) = (9/2)*Pi, (5) = (11/2)*Pi, (6) = (13/2)*Pi, (7) = (15/2)*Pi, (8) = (17/2)*Pi, (9) = (19/2)*Pi, (10) = (21/2)*Pi})

[seq(fsolve(f, q = OddAsymptotes[i] .. OddAsymptotes[i+1]), i = 1 .. n-1)];

[7.853797468, 10.99548650, 14.13711266, 17.27872097, 20.42032239, 23.56192056, 26.70351698, 29.84511237, 32.98670709]

Vector[row](n-1, proc (i) options operator, arrow; fsolve(f, q = OddAsymptotes[i] .. OddAsymptotes[i+1]) end proc)

Vector[row](9, {(1) = 7.853797468, (2) = 10.99548650, (3) = 14.13711266, (4) = 17.27872097, (5) = 20.42032239, (6) = 23.56192056, (7) = 26.70351698, (8) = 29.84511237, (9) = 32.98670709})

 

``

Download Asymptotes_acc.mw

Isn't this a repeat of one of your previous Questions?

For larger (symmetric) examples you will get a result more quickly by computing the eigenvalues of the Matrix (cast suitably to float[8] datatype, symmetric indexing function, etc).

You wrote in a followup comment, "In general, the matrices I encounter are real symmetric", so below I'll focus on code for that case.

Here's code for doing it via Eigenvalues, and an example where the SingularValues computation takes about 3.2 times as much wall-clock time, on my Linux machine. (That slowdown ratio could typically grow, as the Matrix size increases, but it can depend on the nature of the Graph data.)

SR := proc(G::Graph)
  local M := GraphTheory:-AdjacencyMatrix(G);
  max(abs~(LinearAlgebra:-Eigenvalues(
    rtable(`if`([rtable_indfns(M)]=[':-symmetric'],
           ':-symmetric',NULL), M,
           'datatype'=':-hfloat','storage'=':-rectangular',
           'order'=':-Fortran_order'))));
end proc:

M22:= GraphTheory:-SpecialGraphs:-M22Graph();

   M22 := Graph 1: an undirected unweighted graph with
          77 vertices and 616 edge(s)

CodeTools:-Usage(SR(M22), iterations=100);

memory used=270.28KiB, alloc change=38.59MiB, cpu time=1.84ms,
real time=1.86ms, gc time=110.19us

            16.0000000000000

G := GraphTheory:-SpecialGraphs:-CirculantGraph(2345,[1,2]);

    G := Graph 2: an undirected unweighted graph with
         2345 vertices and 4690  edge(s)

CodeTools:-Usage(SR(G));

memory used=242.85MiB, alloc change=107.52MiB, cpu time=5.75s,
real time=2.73s, gc time=312.03ms

            4.00000000000000

SpectralRadius_ac.mw

If you know certain additional qualities of the Graph up front (eg. Bipartite) then you might be able to attain even better speed by computing just the "first" selected eigenvalue. I had an old post on computing only selected eigenvectors/values. But it's not helpful if you have to call it twice to compute first and last (largest-positive and smallest-negative). I'm not sure whether this is workable.

These are all using external calling to LAPACK for dense matrix operations. ie. not sparse stuff.

Are you asking about how to get the exploration to appear in a new document, in a separate GUI tab? If so, then it can be done by using the Explore command explicitly, instead using rigght-click context-menu.

For example,

restart;

Explore(plot(sin(a*x)*x^2 + cos(b*x),
             x = -2*Pi .. 2*Pi, view = -30 .. 30),
        a = 1.0 .. 10.0, b = 1.0 .. 10.0,
        animate, newsheet, showbanner = false);

 

If you mean something else then you really ought to explain it more clearly. What do you mean by "save", and "extension file"?

That particular accent you've chosen may be getting parsed (in 2D Input) as conjugate(a). So you might be seeing a 2D Input version of this:

conjugate(a):=3;

3

conjugate(a)+2;

conjugate(a)+2

Download conjugate_sigh.mw

The same kind of thing seems to happen in the 2D Input parsing if I utilize the Layout palette (for "Over") and the Punctuation palette for inserting the "macr" over-bar. That is, I get a call to conjugate as the parsed meaning.

You could get the expected behaviour using the "Over" item from the Layout palette and then some other kind of over-bar from the Punctuation palette, eg. the "mdash", "minus", or "ndash" items. As over-bars those don't seem to make the names get parsed as calls to conjugate, instead parsing as distinct (typeset) names that behave more as you'd expected.

Please check whether you have shown the differential equation, and the value of lambda, correctly. (For example, what do you expect for 4*f(8) = x*f(eta) when x=4 and eta=8?)

As I parse the DE that you've shown, the value of 4*f(8) is only about 3.2 when lambda=0.5. But that doesn't correspond to the red color for the values near 4.0 in your colorbar. So either I've misunderstood the DE, or something else is not right.

Below I use lambda=0.005.

restart;

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

Typesetting:-Settings(typesetprime=true):
Typesetting:-Suppress(f(x));
interface(typesetting=extended):

de := diff(f(x),x,x,x) + (f(x)*diff(f(x),x,x) - diff(f(x),x)^2)
      - lambda*diff(f(x),x) = 0;

diff(diff(diff(f(x), x), x), x)+f(x)*(diff(diff(f(x), x), x))-(diff(f(x), x))^2-lambda*(diff(f(x), x)) = 0

lambda := 0.005:

desol := dsolve({de, f(0)=0, D(f)(0)=1, D(f)(10)=0}, f(x),
                numeric, output=listprocedure):

fsol := eval(f(x), desol):

4 * fsol(8);

HFloat(3.987786123641983)

plots:-display(
  plots:-densityplot(x*fsol(y), x=0..4, y=0..8,
                     style=surface, restricttoranges,
                     colorscheme=["zcoloring",
                                  z->2/3*(1-z/4)]),
  plots:-contourplot(x*fsol(y), x=0..4, y=0..8,
                     contours=[seq(0..4,0.5)], thickness=0),
  seq(plottools:-line([0,0],[0,0],thickness=10,
                      color=ColorTools:-Color("HSV",[2/3*(4-i)/(4-0),1,1]),
                      legend=sprintf("%.3f",i)),
      i=4..0,-0.5),
  labels=[x,eta], labeldirections=[horizontal,vertical],
  legendstyle=[location=right], axes=box,
  axis[1]=[tickmarks=[seq(0.0..4.0,0.5)]],
  axesfont=[Times,14], size=[500,400]
);

Download de_col.mw

 

[edit] Looking at the OP's previous postings it seems as if he might be using the significantly older Maple 17, in which case some of the above functionality is not available. Here is a revision, for that older version.

My query about the correctness of the DE and value of lambda still applies.

restart;

kerneopts(version);

kerneopts(version)

de := diff(f(x),x,x,x) + (f(x)*diff(f(x),x,x) - diff(f(x),x)^2)
      - lambda*diff(f(x),x) = 0;

diff(diff(diff(f(x), x), x), x)+f(x)*(diff(diff(f(x), x), x))-(diff(f(x), x))^2-lambda*(diff(f(x), x)) = 0

lambda := 0.005:

desol := dsolve({de, f(0)=0, D(f)(0)=1, D(f)(10)=0}, f(x),
                numeric, output=listprocedure):

fsol := eval(f(x), desol):

4 * fsol(8);

HFloat(3.9877861236419827)

plots:-display(
  plots:-densityplot(piecewise( x<-0.75, 0, x<-0.5, 1, x<0, 2/3,
                                2/3*(1 - x*fsol(y)/4) ),
                     x=-1..4, y=0..8,
                     style=patchnogrid,
                     colorstyle=HUE),
  plots:-contourplot(x*fsol(y), x=0..4, y=0..8,
                     contours=[seq(0..4,0.5)], thickness=0),
  seq(plottools:-line([0,0],[0,0],thickness=10,
                      color=ColorTools:-Color("HSV",[2/3*(4-i)/(4-0),1,1]),
                      legend=sprintf("%.3f",i)),
      i=4..0,-0.5),
  labels=[x,eta], labeldirections=[horizontal,vertical],
  legendstyle=[location=right], axes=box, view=[0..4,0..8],
  axis[1]=[tickmarks=[seq(0.0..4.0,0.5)]],
  axesfont=[Times,14]
);

Download de_col17.mw

You are substituting 4 for instances of 3.

subs(3=4,3);

             4

subs(3=4,9);

             9

x:=3;

             3

subs(x=4, 3);

             4

subs(x=4, 9);

             9

Check whether this (or any other Answer) provides as many roots as expected. (Eg, compare against the plot, or compare methods against each other, etc.)

restart:

with(Student:-Calculus1):

a[1] := .1093: k[3] := 7.5*10^(-12): k[2] := 3.8*10^(-12):
d := 0.2e-3: eta[1] := 0.240e-1: alpha[2] := -.1104:
alpha[3] := -0.1104e-2: eta[2] := .1361: xi := 1.219*10^(-6):
alpha := 1-alpha[3]^2/(a[1]*eta[1]):theta[0] := 0.5e-1:
Hc := (Pi/d)*sqrt(k[2]/xi): H := 5.5*Hc: lambda := a[1]/(xi*H^2):

 

sols := Roots((H/(Hc))^2 -(4*q^2)/Pi^2*((tan(q)- q/(1-alpha))/(tan(q)-q)),
              q=0..100, numeric, maxsols=200);

[1.579423766, 4.712622100, 7.853994064, 10.99555811, 14.13714635, 17.27873942, 20.42033357, 23.56192783, 26.70352198, 29.84511595, 32.98670974, 36.12830339, 39.26989691, 42.41149032, 45.55308364, 48.69467689, 51.83627007, 54.97786320, 58.11945627, 61.26104931, 64.40264231, 67.54423529, 70.68582823, 73.82742115, 76.96901405, 80.11060693, 83.25219980, 86.39379265, 89.53538549, 92.67697831, 95.81857112, 98.96016393]

nops(sols);

32

for i from 1 to nops(sols) do
 p[i] := alpha*lambda/(1-(4*sols[i])^2*(Hc/H)^2/Pi^2);
end do;

4.446263879

-20.25154621

-1.670508907

-.7029799724

-.3966600833

-.2567907775

-.1804394088

-.1339680365

-.1035031442

-0.8241813528e-1

-0.6720612510e-1

-0.5586412774e-1

-0.4717829446e-1

-0.4037722719e-1

-0.3495120384e-1

-0.3055223066e-1

-0.2693603757e-1

-0.2392701066e-1

-0.2139623682e-1

-0.1924733487e-1

-0.1740706641e-1

-0.1581896595e-1

-0.1443892890e-1

-0.1323209774e-1

-0.1217062731e-1

-0.1123205608e-1

-0.1039810281e-1

-0.9653765777e-2

-0.8986640837e-2

-0.8386399464e-2

-0.7844385207e-2

-0.7353298899e-2

 

Download rf_examp.mw

First 71 72 73 74 75 76 77 Last Page 73 of 336