pagan

5147 Reputation

23 Badges

17 years, 127 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are answers submitted by pagan

It looks like fsolve was getting confused about the dummy variable of integration, _t0, in that symbolic CDF result. That looks like a bug.

There is more than one workaround, I think. And those can differ in their respective speeds of execution and performance.

> with(Statistics):

> X := RandomVariable(NonCentralFRatio(nu, omega, delta)):
> nu:=3;omega:=3*n; delta:=4/(2*100)*(2.5^2+1.5^2+3.5^2+4.5^2):

> F:=N->CDF(RandomVariable(NonCentralFRatio(op(eval(
>              [nu,omega,delta],n=N)))),100,numeric)-0.9:

> st:=time(): fsolve(F,0..1); time()-st;
                                 0.3031007062

                                     0.924

> # plot(F,0..1);

Compare these timings, for just a single given value of n,

> st:=time():
> evalf(eval(CDF(X,100),n=0.67))-0.9;
                                 0.0876176466

> time()-st;
                                     2.212

> st:=time():
> F(0.67);
                                 0.0876176469

> time()-st;
                                     0.004

B is a Vector, not a Matrix, for that to make sense, yes?

> restart:
> dat:=DEtools[convertsys]({diff(y(x),x$3)-2*diff(y(x),x$2)
>                           +x^2*diff(y(x),x)+sin(x)*y(x)-1=0},
>                          {},y,x,'yv','ypv'):

> XP:=Vector(3,symbol=ypv):
> X:=Vector(3,symbol=yv):

> A,XPB:=LinearAlgebra:-GenerateMatrix(map(rhs-lhs,dat[1]),
>                                      convert(X,list)):

> B:=XP-XPB:

> XP=A.X+B;
              [ypv[1]]   [                yv[2]                ]
              [      ]   [                                     ]
              [ypv[2]] = [                yv[3]                ]
              [      ]   [                                     ]
              [ypv[3]]   [               2                     ]
                         [1 + 2 yv[3] - x  yv[2] - sin(x) yv[1]]
 
> dat[1];
[ypv[1] = yv[2], ypv[2] = yv[3],
 
                            2
    ypv[3] = 1 + 2 yv[3] - x  yv[2] - sin(x) yv[1]]

> XP,A,X,B;
                 [ypv[1]]  [   0        1     0]  [yv[1]]  [0]
                 [      ]  [                   ]  [     ]  [ ]
                 [ypv[2]], [   0        0     1], [yv[2]], [0]
                 [      ]  [                   ]  [     ]  [ ]
                 [ypv[3]]  [             2     ]  [yv[3]]  [1]
                           [-sin(x)    -x     2]

How is the data supplied by the website?

Joe and Markiyan have already shown you how to get the result most simply, using a stock routine.

But if you are trying to teach yourself Maple, you might still be interested in programming it yourself. Now, usually you can print the contents of Maple Library routines. But `taylor` is a (compiled) kernel builtin and not so visible.

Here are two simple (read: crude) approaches. The first is for f an expression and the second is for f a procedure (of one argument).

> Taylor:=proc(f,t,a,n)
> local i;
>   eval(f,t=a)+add(eval(diff(f,t$i),t=a)/(i!)*(t-a)^(i),i=1..n);
> end proc:
> Taylor(ln(x),x,1,4);
                                   2          3          4
                            (x - 1)    (x - 1)    (x - 1)
                    x - 1 - -------- + -------- - --------
                               2          3          4
 
> Taylor:=proc(f,t,a,n)
> local i;
>   unapply(add((D@@i)(f)(a)/(i!)*(t-a)^(i),i=0..n),t);
> end proc:
> Taylor(ln,x,1,4);
                                    2              3              4
            x -> x - 1 - 1/2 (x - 1)  + 1/3 (x - 1)  - 1/4 (x - 1)

Of course, you could also use unapply on the first result, to get the (operator) second result.

Note that `diff` has `option remember` and thus a remember table to store previously computed results. This can make it more efficient when computing lots of derivatives (of various orders) of an expression.

For example

> other_eqs:={y=cos(t)-5.0, z=tan(t)+9.1}:

> x:=[0.1, 0.2, 0.3]:

> f:=sin(t)-X^2:

> for i in x do
>   s(i):=fsolve({eval(f,X=i)},{t});
> end do;
                         s(0.1) := {t = 0.01000016667}
 
                         s(0.2) := {t = 0.04001067435}
 
                         s(0.3) := {t = 0.09012194501}
 
> eval(other_eqs, s(0.2)); # or s(0.1), etc
                      {y = -4.000800320, z = 9.140032038}
 
> for i from 1 to nops(x) do
>   s(i):=fsolve({eval(f,X=x[i])},{t});
> end do;
                          s(1) := {t = 0.01000016667}
 
                          s(2) := {t = 0.04001067435}
 
                          s(3) := {t = 0.09012194501}
 
> eval(other_eqs, s(2)); # or s(1), etc
                      {y = -4.000800320, z = 9.140032038}

The methodology would look the same, even if you had more than one unknown t in the fsolve calls (and hence also in the solutions s(i) which are returned and assigned as sets).

So, there are two ways above. One has i take values from x, while the other has i take ordinal values to denote the position in x.

BTW, I used f as an expression. But you might have used f as an operator, passing f(i) of f(x[i]) respectively to fsolve.

If you are looking for a floating-point approximation to the value of t, then look at the fsolve command.

If you are looking for an exact result, then you might look at the solve command.

You wrote "function", which might mean (to you) either a Maple expression or a Maple procedure. Below is fsolve working on each of those.

> restart: # first, an expression u

> u := sin(t)+t^2;
                                              2
                               u := sin(t) + t
 
> fsolve(u=1.23456789, t=0..3);
                                 0.7457059272

> eval(u,t=%);
                                  1.234567890
 
> restart: # next, an operator u

> u := t -> sin(t)+t^2;
                                                 2
                             u := t -> sin(t) + t
 
> fsolve(u=1.23456789, 0..3);
                                 0.7457059272
 
> fsolve(u(t)=1.23456789, t=0..3);
                                 0.7457059272

> u(%);
                                  1.234567890

Remove the space between the 5 and the f in your `%12.5 f`.

Presumably this is just a simplified example of your actual task, because for something this short and simple it would be easier to use expressions instead of procedures.

aold:=x-> x^2;

for i from 1 to 5 do
   anew:=unapply(aold(x)+x,x);
   aold:=eval(anew);
end do;

aold:='aold':
eval(aold);
eval(anew); # there is no dependency upon aold

You might be better off using expressions instead of procedures, regardless of how complicated the actual task is. With expressions, it's much easier.

restart:

aold:=x^2;
 
for i from 1 to 5 do
   anew:=aold+x;
   aold:=anew;
end do;

aold:='aold':
aold;
anew;

eval(anew,x=17.3); # you don't need a function, to get a value!

f_anew:=unapply(anew,x): # ...but maybe you really prefer procedures
f_anew(17.3);
Try the help page ?fnormal That should let you filter by floating-point magnitude (as opposed to by imaginariness/realness).

A table is of type last_name_eval (see also here).

You can get the result 'table' by querying whattype(eval(T,1)) for table T.

This doesn't produce exactly what you described. There isn't a stage with | 4 - 9 |

Also, it takes a dubious shortcut. Instead of reliably working from the inside out, in terms of depth, it just presumes that the ordering of subexpressions by depth matches up with the ordering by length.

A more robust way might be to treat, at any given stage, all (type specfunc) _Inert_PROD, _Inert_SUM,  _Inert_POWER, or _Inert_RATIONAL function calls whose arguments are calls of _Inert_NAME, _Inert_INTPOS, or _Inert_INTNEG.

Maybe someone will show an even neater way.

> restart:

> F:=X->subsindets(X,identical(FromInert(
      sort(convert(indets(ToInert(X),
      specfunc(anything,{_Inert_PROD,_Inert_SUM})),list),
      length)[1])),t->subsindets(t,name,tt->parse(tt))):

> expr:=`12`/(`1`-abs(`2`^2-`3`^2));
                                          12
                            expr := ---------------
                                           2    2
                                    1 - | 2  - 3  |
 
> expr:=F(expr);
                                          12
                            expr := --------------
                                           2
                                    1 - | 2  - 9 |
 
> expr:=F(expr);
                                          12
                              expr := ----------
                                      1 - | -5 |
 
> expr:=F(expr);
                                           12
                                expr := - ----
                                           4
 
> expr:=F(expr);
                                  expr := -3

Your second version does not place any assumptions sigma.

> restart:

> assume(0 < sigma);

> is(sigma>0);
                                     true
 
> about(sigma);
Originally sigma, renamed sigma~:
  is assumed to be: RealRange(Open(0),infinity)
 
> restart:

> 0 < sigma;
                                   0 < sigma
 
> is(sigma>0);
                                     false
 
> about(sigma);
sigma:
  nothing known about this object

Perhaps you were thinking of something like this.

> restart:

> ineq := 0 < sigma;
                               ineq := 0 < sigma
 
> is(sigma>0) assuming ineq;
                                     true

In the first method, try replacing

   X2:= min(map(rhs,R));

   X1:= max(map(rhs,R));   

with

   X2:= min(op(map(rhs,R)));

   X1:= max(op(map(rhs,R)));   

It seems due to this older behaviour of Maple 10

> min({1.2,2.3});
Error, (in simpl/min) arguments must be of type algebraic

> min(op({1.2,2.3}));
                                      1.2

while in Maple 13,

> min({1.2,2.3});
                                      1.2

Suppose that you have read your data into a 2xN Matrix (where N may vary). You can create a `piecewise` dynamically from that Matrix as follows.

> A:=LinearAlgebra:-RandomMatrix(2,3);
                                 [57    -93    -72]
                            A := [                ]
                                 [27    -76     -2]
 
> N:=op([1,2],A);
                                    N := 3
 
> f:=piecewise(seq(op([x<A[1,j],A[2,j]]),j=1..op([1,2],A)),666);
                               { 27          x < 57
                               {
                               { -76         x < -93
                          f := {
                               { -2          x < -72
                               {
                               { 666        otherwise
 
> eval(f,x=-100);
                                      27

Note, however, that evaluation of `f` happens in a particular first-satisfied-first-returned kind of way. If the conditions of the piecewise are not "sorted" then you may not get what you hoped for. Using the above `f` you may not be able to attain the -76 and -2 result.

> eval(f,x=-80);
                                      27

So you might wish to sort the piecewise data. (The following is not even close to most efficient, for large amounts of data. Perhaps search for mapleprimes threads on sorting with attributes.)

> As:=Matrix(sort([seq(A[1..-1,j],j=1..op([1,2],A))],
>                 (a,b)->`if`(a[1]<b[1],true,false)));
                                 [-93    -72    57]
                           As := [                ]
                                 [-76     -2    27]
 
> fs:=piecewise(seq(op([x<As[1,j],As[2,j]]),j=1..op([1,2],As)),666);
                               { -76         x < -93
                               {
                               { -2          x < -72
                         fs := {
                               { 27          x < 57
                               {
                               { 666        otherwise
 
> eval(fs,x=-80);
                                      -2

You didn't mention whether all your data was numeric (as opposed to symbolic), or how much of it there was, or whether your conditionals (breakpoints) would be as simple univariate inequality tests, etc, etc. If you have a LOT of numeric data with such simple conditionals then you might get much more efficient evaluation (...at many x points, further down the road in your calculations) if you used ArrayInterpolation or rolled your own (possibly Compile'able) procedure instead of using `piecewise`.

Here's the code to my original reply (the one that got accidentally replaced by Thomas's reply entitled "kernelopts(datadir)".

It's not fast, but I'm sure it could be made faster, and smarter. The basic idea is that an image, once read into Maple, can be put up as either a 3D or a 2D plot.

First the 3D plot. (Thanks, Thomas, for the datadir tip.)

restart:
a,b:=-10,10:
P1:=plots:-spacecurve([sin(t),t,10],t=a..b,color=black):
#immy:=ImageTools:-Read("/usr/local/maple10.03/data/images/phone.jpg"):
immy:=ImageTools:-Read(cat(kernelopts(datadir),"/images/phone.jpg")):
m,n:=op(map2(op,2,[rtable_dims(immy)])):
P2a:=ImageTools:-Preview(immy):
P2:=plottools:-translate(P2a,-floor(m/2),-floor(n/2),1):
P2:=plottools:-scale(P2,(b-a)/m,(b-a)/n,1):
plots:-display([P2,P1],orientation=[0.0,0.1]);

Next, the 2D plot, using P2a from above.

Z:=op(2,op(4,op(1,P2a))):
Zc:=Array(1..265*265,1..3,(i,j)->Z[265-max(1,iquo(i,265))+1,265-irem(i,265),j],datatype=float[8]):
M:=Matrix(265*265,2,datatype=float[8]):
V:=Vector(265,i->265-i+1,datatype=float[8]):
Vc:=Vector(265,datatype=float[8]):
for i from 1 to 265 do
  ArrayTools:-Fill(265,i,Vc,0,1);
  M[(i-1)*265+1..(i-1)*265+265,2]:=Vc:
  M[(i-1)*265+1..(i-1)*265+265,1]:=V:
end do:
P3a:=plots:-pointplot(convert(M,listlist),color=[seq(COLOR(RGB,seq(Zc[i][k],k=1..3)),i=1..265^2)]):
P3:=plottools:-translate(P3a,-floor(m/2),-floor(n/2)):
P3:=plottools:-scale(P3,(b-a)/m,(b-a)/n):
P4:=plot(sin(x),x=a..b,color=black):
plots:-display([P3,P4]);

Anyway, it can be done, is the essence.

First 34 35 36 37 38 39 40 Last Page 36 of 48