acer

32373 Reputation

29 Badges

19 years, 333 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Three popular choices for trying to find all roots in a finite real range are to use Student:-Calculus1:-Roots (which calls fsolve repeatedly with its `avoid` option), RootFinding:-Analytic, or to repeatedly call RootFinding:-NextZero.

Below, the original problem is used, but for a larger domain of x=0 to x=100.

Note that, with default settings at least, the fsolve approach finds fewer solutions. I would guess that for a large enough right end-point the extreme steepness of the mathematical function between the upper roots would cause `NextZero` to falter as well. It may not be an issue for this particular problem, but the RootFinding:-Analytic approach can messy if there are roots very close to the real line, since it can be difficult to distinguish which are valid but unwanted nonreal solutions and which are wanted purely real solutions which just happen to have been approximated with very small but nonzero imaginary components due to roundoff error.

restart:

findroots:=proc(expr,a,b,{guard::posint:=5,maxtries::posint:=50})
local F,x,sols,i,res,start,t;
   x:=indets(expr,name) minus {constants};
   if nops(x)>1 then error "too many indeterminates"; end if;
   F:=subs(__F=unapply(expr,x[1]),__G=guard,proc(t)
      Digits:=Digits+__G;
      __F(t);
   end proc);
   sols,i,start:=table([]),0,a;
   to maxtries do
      i:=i+1;
      res:=RootFinding:-NextZero(F,start,
                                 'maxdistance'=b-start);
      if type(res,numeric) then
         sols[i]:=fnormal(res);
         if sols[i]=sols[i-1] then
            start:=sols[i]+1.0*10^(-Digits);
            i:=i-1;
         else
            start:=sols[i];
         end if;
      else
         break;
      end if;
   end do;
   op({entries(sols,'nolist')});
end proc:

CodeTools:-Usage( findroots(exp(x)*cos(x)+1, 0, 100) );
memory used=0.55MiB, alloc change=0 bytes, cpu time=47.00ms, real time=33.00ms

1.746139530, 4.703323759, 7.854369686, 10.99555751, 14.13716766, 17.27875956, 

  20.42035224, 23.56194490, 26.70353755, 29.84513020, 32.98672286, 

  36.12831551, 39.26990816, 42.41150082, 45.55309347, 48.69468613, 

  51.83627878, 54.97787143, 58.11946409, 61.26105674, 64.40264939, 

  67.54424205, 70.68583470, 73.82742735, 76.96902001, 80.11061266, 

  83.25220532, 86.39379797, 89.53539062, 92.67698328, 95.81857593, 98.96016858

nops({%});
                                     32
restart:

CodeTools:-Usage( RootFinding[Analytic](exp(x)*cos(x)+1=0, re=0..100, im=-1..1) );
memory used=23.34MiB, alloc change=24.00MiB, cpu time=343.00ms, real time=341.00ms

 48.6946861306418, 45.5530934770530, 42.4115008234622, 39.2699081698730, 

   29.8451302091029, 73.8274273593600, 70.6858347057710, 67.5442420521805, 

   64.4026493985910, 61.2610567450010, 86.3937979737195, 83.2522053201305, 

   80.1106126665395, 76.9690200129500, 92.6769832808990, 89.5353906273097, 

   95.8185759344885, 98.9601685880785, 54.9778714378215, 51.8362787842320, 

   58.1194640914110, 23.5619449018649, 10.9955575115013, 20.4203522496875, 

   32.9867228626928, 26.7035375555158, 36.1283155162826, 7.85436968657417, 

   4.70332375945224, 14.1371676661008, 17.2787595634161, 1.74613953040801


nops([%]);
                                     32

restart:

asolve := proc(ex, var :: name, rng :: range, {avoid :: set := {}})
local avoids,ends,i,sols;
    sols := [fsolve](ex, var, rng, _options['avoid']);
    if sols = [] or not sols :: 'list(numeric)' then
        return avoid;
    else
        ends := [lhs(rng), op(sort(sols)), rhs(rng)]; # sols already sorted?
        avoids := {op(avoid),map2(`=`,var,sols)[]};
        {seq(op(thisproc(ex,var,ends[i-1]..ends[i], ':-avoid' = avoids))
             , i = 2 .. nops(ends))};
    end if;
end proc:

CodeTools:-Usage( asolve(exp(x)*cos(x)+1, x, 0..100) );
memory used=89.20MiB, alloc change=28.00MiB, cpu time=1.36s, real time=1.37s

   {x = 1.746139530, x = 4.703323759, x = 7.854369687, x = 10.99555751, 

     x = 14.13716767, x = 17.27875956, x = 20.42035225, x = 23.56194490, 

     x = 26.70353756, x = 29.84513021, x = 32.98672286, x = 36.12831552, 

     x = 39.26990817, x = 42.41150082, x = 45.55309348, x = 48.69468613, 

     x = 51.83627878}

nops(%);
                                     17
restart:
CodeTools:-Usage( Student:-Calculus1:-Roots(exp(x)*cos(x)+1,x=0..100,numeric) );
memory used=87.17MiB, alloc change=28.00MiB, cpu time=1.36s, real time=1.36s

[1.746139530, 4.703323759, 7.854369687, 10.99555751, 14.13716767, 17.27875956, 

  20.42035225, 23.56194490, 26.70353756, 29.84513021, 32.98672286, 

  36.12831552, 39.26990817, 42.41150082, 45.55309348, 48.69468613, 51.83627878

  ]

nops(%);
                                     17

acer

You could try using the Basis command instead.

with(LinearAlgebra):

A:=Matrix([[-3,6,-1,1-7],[1,-2,2,3,-1],[2,-4,5,8,-4]]);

                               [-3   6  -1  -6   0]
                               [                  ]
                          A := [ 1  -2   2   3  -1]
                               [                  ]
                               [ 2  -4   5   8  -4]

Basis({seq(A[..,i],i=1..ColumnDimension(A))});

                             /[-3]  [-1]  [-6]\ 
                             |[  ]  [  ]  [  ]| 
                            < [ 1], [ 2], [ 3] >
                             |[  ]  [  ]  [  ]| 
                             \[ 2]  [ 5]  [ 8]/ 

Use the menu item Edit -> Split or Join.

It has the keyboard shortcut of F3 (on MS-Windows at least).

acer

Did you intend y(x)^M instead of y^M?

restart:
eq1:=diff(y(x),x$2)+2/x*(diff(y(x),x))+y(x)^M=0:
ic:=y(0)=a,D(y)(0)=0:
p:=dsolve(eval({eq1,ic},[a=1,M=3]), y(x),numeric,output=listprocedure):
Y:=eval(y(x),p):
plot(Y,0..500);

Two basic approaches to plotting a surface from 3D data are interpolation and smoothing. In the former, the surface goes through the data points exactly.

Your "grid" of data points is not full. There are gaps. If you had a full grid (lattice) of x-z data points then plots:-surfdata would handle it quickly, to interpolate a surface.

Stock Maple does not presently have a command to interpolate (higher than linearly) an irregularly spaced set of 3D data points to get a smooth surface. The `surfdata` command can produce a collection of polygon plots as a linear interpolation. But it's not very attractive, in part because of shading and partly because it's only piecewise smooth.

I converted your data to csv, after chopping off the first row of the spreadsheet.

restart:
M:=ImportMatrix("y0.csv",source=csv,datatype=float[8]):

For reasons I don't fully understand, `surfdata` sometimes stalls while trying to handlt the full 383x3 Matrix. I split it into parts, manually. (.. and even then, sometimes it would hang, and then immediately retrying the problematic call would immediately succeed. Weird.)

P1:=plots:-surfdata(M[1..162],0.0..0.4,0.05..1.0):
P2:=plots:-surfdata(M[142..225],0.45..0.55,0.0..1.0):
P3:=plots:-surfdata(M[205..315],0.60..0.80,0.0..1.0):
P3a:=plots:-surfdata(M[298..336],0.80..0.85,0.0..1.0):
P4:=plots:-surfdata(M[316..],0.85..1.0,0.0..1.0):
plots:-display(P1,P2,P3,P3a,P4);

A simpler, if slower and more computationally intensive approach is to do "lowess" smoothing. This has quite a few options to it, but the following seems not too bad, and the outliers don't seem to be painfully far from the surface.

Statistics:-ScatterPlot3D(M,lowess=true,bandwidth=0.15,fitorder=2,showpoints=true);

Yet another approach would be interpolate repeatedly in just 1 dimension, using the original data, so as to form a full grid (Matrix). For example, for just one portion, to turn this,

M[110..125,1..3];
         [0.300000000000000                  0.  -60.3200000000000]
         [                                                        ]
         [0.300000000000000  0.0500000000000000  -62.8000000000000]
         [                                                        ]
         [0.300000000000000   0.100000000000000  -63.6600000000000]
         [                                                        ]
         [0.300000000000000   0.150000000000000  -71.2700000000000]
         [                                                        ]
         [0.300000000000000   0.200000000000000  -60.2500000000000]
         [                                                        ]
         [0.300000000000000   0.250000000000000  -56.1800000000000]
         [                                                        ]
         [0.300000000000000   0.300000000000000  -49.9900000000000]
         [                                                        ]
         [0.300000000000000   0.500000000000000  -53.5200000000000]
         [                                                        ]
         [0.300000000000000   0.550000000000000  -47.0100000000000]
         [                                                        ]
         [0.300000000000000   0.600000000000000  -57.0500000000000]
         [                                                        ]
         [0.300000000000000   0.650000000000000  -56.8900000000000]
         [                                                        ]
         [0.300000000000000   0.700000000000000  -56.7300000000000]
         [                                                        ]
         [0.300000000000000   0.750000000000000  -55.5200000000000]
         [                                                        ]
         [0.300000000000000   0.850000000000000  -50.2100000000000]
         [                                                        ]
         [0.300000000000000   0.950000000000000  -56.4100000000000]
         [                                                        ]
         [0.300000000000000                  1.  -60.3200000000000]

into this,

Matrix([Vector(21,0.3),
        Vector(21,i->(i-1)*0.05),
        CurveFitting:-ArrayInterpolation(M[110..125,2..3],
                                         Vector(21,i->(i-1)*0.05))],
        datatype=float[8]);
         [0.300000000000000                  0.  -60.3200000000000]
         [                                                        ]
         [0.300000000000000  0.0500000000000000  -62.8000000000000]
         [                                                        ]
         [0.300000000000000   0.100000000000000  -63.6600000000000]
         [                                                        ]
         [0.300000000000000   0.150000000000000  -71.2700000000000]
         [                                                        ]
         [0.300000000000000   0.200000000000000  -60.2500000000000]
         [                                                        ]
         [0.300000000000000   0.250000000000000  -56.1800000000000]
         [                                                        ]
         [0.300000000000000   0.300000000000000  -49.9900000000000]
         [                                                        ]
         [0.300000000000000   0.350000000000000  -50.8725000000000]
         [                                                        ]
         [0.300000000000000   0.400000000000000  -51.7550000000000]
         [                                                        ]
         [0.300000000000000   0.450000000000000  -52.6375000000000]
         [                                                        ]
         [0.300000000000000   0.500000000000000  -53.5200000000000]
         [                                                        ]
         [0.300000000000000   0.550000000000000  -47.0100000000000]
         [                                                        ]
         [0.300000000000000   0.600000000000000  -57.0500000000000]
         [                                                        ]
         [0.300000000000000   0.650000000000000  -56.8900000000000]
         [                                                        ]
         [0.300000000000000   0.700000000000000  -56.7300000000000]
         [                                                        ]
         [0.300000000000000   0.750000000000000  -55.5200000000000]
         [                                                        ]
         [0.300000000000000   0.800000000000000  -52.8650000000000]
         [                                                        ]
         [0.300000000000000   0.850000000000000  -50.2100000000000]
         [                                                        ]
         [0.300000000000000   0.900000000000000  -53.3100000000000]
         [                                                        ]
         [0.300000000000000   0.950000000000000  -56.4100000000000]
         [                                                        ]
         [0.300000000000000                  1.  -60.3200000000000]

If that were done for all the sets of x-values, then (only) all the final y-column results could be put into just a single Matrix which could be fed to `surfdata` (and the x- and z-ranges supplied as simple range arguments). I haven't done this.

acer

There are several other ways. Here are two.

a:={1,2}:

Vector([a[]]);
                                     [1]
                                     [ ]
                                     [2]

<a[]>;
                                     [1]
                                     [ ]
                                     [2]

But note that you may not always get the ordering of entries that you expect, which in turn could possibly make coding the construction of your Matrices difficult.

restart:

a:={2,11,-1}:
<a[]>;
                                    [-1]
                                    [  ]
                                    [ 2]
                                    [  ]
                                    [11]

a:={x,c,X}:
<a[]>;
                                     [X]
                                     [ ]
                                     [c]
                                     [ ]
                                     [x]

Some routines can accept certain of their arguments in either `list` or `set` form, where the ordering of names may be repected when using `list` input. I am not sure whether that fact would help your case. An example is `solve`,

restart:

solve({a+b-4,a-b-1},[a,b]);
                              [[    5      3]]
                              [[a = -, b = -]]
                              [[    2      2]]

solve({a+b-4,a-b-1},[b,a]);
                              [[    3      5]]
                              [[b = -, a = -]]
                              [[    2      2]]

Sometimes there are other ways around that, including 2-argument eval.

restart:

sol := solve({a+b-4,a-b-1},{a,b});
                                   /    5      3\ 
                           sol := { a = -, b = - }
                                   \    2      2/ 

eval([b,a],sol);
                                   [3  5]
                                   [-, -]
                                   [2  2]

It might be easier to help more if you could show a more detailed problematic example.

acer

> restart:
> d:=evalf(Pi,1);         
                                    d := 3.

> evalf(subs(x=d,1/x),1); 
                                      0.3

> evalf(subs(x=d,3/x),1);
                                      0.9

> restart:
> d:=evalf(Pi,2);        
                                   d := 3.1

> evalf(subs(x=d,1/x),2);
                                     0.32

> evalf(subs(x=d,3/x),2);
                                     0.96

> restart:               
> d:=evalf(Pi,3);        
                                   d := 3.14

> evalf(subs(x=d,1/x),3);
                                     0.318

> evalf(subs(x=d,3/x),3);
                                     0.954

And, for interest, following this,

 restart:
> ee:=y/x:
> d:=evalf(Pi,1);
                                    d := 3.

> evalf(eval(ee,[x=d,y=3]),1);  
                                      1.

> evalf(subs(x=d,y=3,ee),1);
                                      0.9

> restart:
> invee:=x/y:
> ee:=1/invee;
                                   ee := y/x

> d:=evalf(Pi,1);
                                    d := 3.

> evalf(eval(ee,[x=d,y=3]),1);
                                      0.9

> evalf(subs(x=d,y=3,ee),1);
                                      0.9

acer

Did you intend a multiplication sign between the first instance of `lambda` and the opening bracket that follows it?

acer

The command that does what you've described (in your followup note to Carl) is `coeff`.

X:=X_FS2(s)=xp0*s+x0:

XP:=diff(X_FS2(s),s)=xp0:

coeff(op(2,X),x0);
                               1

coeff(op(2,XP),x0);
                               0

acer

My 64bit Maple 17.02 Standard GUI on Windows 7 Pro prints only the empty character box as output when I enter the Maple 1D Notation code,

`&#9813;`;

Can this be reduced to a Maple font issue then?

N:=sscanf(`2655`,"%x")[];
                              9813

convert(N,hex);;
                              2655

Hence I actually tried it as,

cat(`&#`,N,`;`);

and of course that could be sprintf'd instead, to get a string rather than a name.

I see that I can get a text plot involving the square-root symbol if I instead do, say,

plots:-display(plots:-textplot([1,1,"&#8730;"]));

acer

You have assigned an expression to the name `w`, not a procedure. So it's not appropropriate to then use `w` in a function call such as `w(x)`.

Like `plot` and `evalf(Int(...))` and several other commands, `fsolve` can work with either operator form (procedures, say) or expressions. Mixing these up is a common mistake (so don't feel bad about it).

Hence you could be calling fsolve like,

   fsolve(w = 0, x = kk .. kk+1);

instead of as your original,

   fsolve(w(x) = 0, x = kk .. kk+1);

You might get better performance here by using RootFinding:-NextZero instead of fsolve for this problem. For that you would need to produce an operator (procedure) by unapplying expression `w` with respect to name `x`. And then you could apply `NextZero` repeatedly, such as,

wf:=unapply(w,x):
RootFinding:-NextZero(wf,5);
                          6.584620042

RootFinding:-NextZero(wf,%);
                          12.72324078

RootFinding:-NextZero(wf,%);
                          18.95497141

Or you could use it in a fancier way, such as in my answer here.

Note also that your code using `fsolve` may be doing inadvertant duplication of attempts (when it fails to find a root), due to your particular check of whether `zz` is of type `float`. If you are running this code at the top level then your type check call type(zz,float) is going to repeat each failed attempt where `zz` is the unevaluated function call (to itself) that `fsolve` returns. This code would run faster if the check were like type(eval(zz,1), float) or if the computation were run within a procedure having local variable `zz` instead. The point is that full evaluation of the argument `zz` as the argument of the call to `type` will reevaluate those returned calls to fsolve and so also repeat the failing computations. See here for some notes on this.

acer

Change the instance of Textfield[TFxv] to TextField[TFxv] to correct the typo.

acer

Try it as,

int(diff(w(s),s,s),s=0..L,continuous);

You may wish to convert a result involving both `diff` and `D` to just one or the other.

acer

restart:

assume(x>0,1>x):

L:=[x^3,x^4,x^2,x,1,0];

                            [  3    4    2          ]
                       L := [x~ , x~ , x~ , x~, 1, 0]

sort( L ); # the usual, and not what is wanted

                          [            2    3    4]
                          [0, 1, x~, x~ , x~ , x~ ]

sort( L, (a,b)->is(a<b) );

                          [     4    3    2       ]
                          [0, x~ , x~ , x~ , x~, 1]

sort( L, (a,b)->is(a>b) );

                          [         2    3    4   ]
                          [1, x~, x~ , x~ , x~ , 0]

acer

That should be sin(x) not sinx.

acer

First 246 247 248 249 250 251 252 Last Page 248 of 336