acer

9 years, 232 days


These are answers submitted by acer

chebpade

Yesterday at 9:23 AM acer 11072

I don't know whether I accurately transcribed your "solution" given in terms of abbreviated Chebyshev polynomials (which I assign to `OP` below). But here is how it can be converted to rational form, both for your "solution" OP as well as for a better approximation.

I also show the continued fraction form (though I use only the explicit rational polynomial form for the plotting).

restart:

f := cos(x):

fa := numapprox:-chebpade(f, x=0..2, [2,2]);

(HFloat(0.36861925958183306)*T(0, x-1)-HFloat(0.7273082240925188)*T(1, x-1)-HFloat(0.13338631081833283)*T(2, x-1))/(T(0, x-1)+HFloat(0.10915303237358412)*T(1, x-1)+HFloat(0.07088232684241387)*T(2, x-1))

sol := evalf( eval(fa, T='orthopoly[T]') );

(HFloat(1.2293137944926846)-HFloat(0.7273082240925188)*x-HFloat(0.26677262163666565)*(x-1.)^2)/(HFloat(0.819964640784002)+HFloat(0.10915303237358412)*x+HFloat(0.14176465368482774)*(x-1.)^2)

numapprox:-confracform(sol);

-HFloat(1.8817992687356087)-HFloat(3.681482753465016)/(x+HFloat(4.081897985226272)+HFloat(28.466776346835747)/(x-HFloat(5.311938552205301)))

# This OP appears to be what you posted.
OP := (-.221091073962959*T(1, x-1)+.7710737338*T(0, x-1)
       -0.4212446689e-1*T(2, x-1))/
      (0.836360586596837e-1*T(1, x-1)+T(0, x-1)+0.3360079945e-1*T(2, x-1)):

solOP := eval(OP, T='orthopoly[T]' );

(-.221091073962959*x+1.034289275-0.8424893378e-1*(x-1)^2)/(0.836360586596837e-1*x+.8827631418+0.6720159890e-1*(x-1)^2)

numapprox:-confracform(solOP);

-1.253674543-1.729701053/(x+17.66344080+339.4769497/(x-18.41888620))

plot( [cos(x), solOP, sol], x=0..2,
      style=[line,point,point], symbolsize=15,
      symbol=[default, soliddiamond, solidbox],
      adaptive=false, numpoints=20, gridlines=false,
      legend=["cos(x)", "solOP", "sol"] );

 


Download chebpade.mw

acer

D

May 20 2015 acer 11072

The conversion fore and aft (to diff, and then back to D) like Carl suggests seems like a sensible way to proceeed, since as you say you know it works ok for diff form.

I suppose that one could write a program which figured out which modifications of the equation (eq) could be utilized as below. The evaluation at a point (eg, at (x,t)) might also need temporary handling. It seems like more work than fun...

eq := D[2](u) = a^2*D[1,1](u);

                                         2           
                        eq := D[2](u) = a  D[1, 1](u)

expr := D[2,2](u);

                             expr := D[2, 2](u)

eval( expr, [D[2](eq)] ) assuming a::constant;

                               2              
                              a  D[1, 1, 2](u)

eval( %, [D[1,1](eq)] ) assuming a::constant;

                              4                 
                             a  D[1, 1, 1, 1](u)

acer

builtin

May 19 2015 acer 11072

You can use timelimit on a call to int, which is not builtin.

The sentence you've quoted from the timelimit help page is about the fact that the kernel will not poll the time resources whilst inside a call to a kernel builtin. But once control returns from such then the resource query may be done.

The Library command int makes use of kernel builtins, here and there, but this doe not bar you from managing it with timelimit in general.

  restart:

  st,str := time(),time[real]():

  timelimit(3, int((a+a*sin(f*x+e))^(1/2)*(c+d*sin(f*x+e))^(5/2),x) );
Error, (in gcd/DIVIDE) time expired

  time()-st, time[real]()-str;

                          3.136, 3.169

Sure, you might get unlucky and still get the kernel locked in some huge (builtin) expand call, etc, resulting in the kernel not stopping until a longer time interval than your supplied limit has elapsed.

But the fact that that a Library command ever invokes a builtin doesn't by itself imply that you cannot use timelimit to manage it. Also, the kernel will take into account the resources spent inside (returned) builtin calls -- it just won't poll and stop within them. That's my understanding, at least.

acer

This is without "extra parameters". (Sorry to belabor this point, but I'm not sure still that we're talking about the same thing. By "extra parameters" I'm imagining that some of the numeric coefficients in your f[i] expressions are parameters which could vary; in the dsolve,parameters sense.)

Note that the compiled Newtonc is still using J3,f3,s3c as globals, which is not exactly the same as compiling all four procedures together. This is just something to start with...


restart:

Digits:=15;

15

N:=4;

4

f:=array([seq(evalf(-2*x[i]+x[i]^2+0.99/i^2),i=1..N)]);

f := Vector[row](4, {(1) = -2.*x[1]+x[1]^2+.99, (2) = -2.*x[2]+x[2]^2+.247500000000000, (3) = -2.*x[3]+x[3]^2+.110000000000000, (4) = -2.*x[4]+x[4]^2+0.618750000000000e-1})

fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});

{x[1] = .900000000000000, x[2] = .132532421355126, x[3] = 1.94339811320566, x[4] = 0.314314686094742e-1}

eqsA := [seq(F[i]=f[i],i=1..N)]:

irform2 := StatSeq(seq(Assign(op(i)),i=eqsA)):

prccons:= codegen[intrep2maple](Proc(Parameters(x::Array,F::Array),irform2)):

f3:=Compiler:-Compile(prccons):

eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];

[Jac[1, 1] = -2.+2.*x[1], Jac[1, 2] = 0., Jac[1, 3] = 0., Jac[1, 4] = 0., Jac[2, 1] = 0., Jac[2, 2] = -2.+2.*x[2], Jac[2, 3] = 0., Jac[2, 4] = 0., Jac[3, 1] = 0., Jac[3, 2] = 0., Jac[3, 3] = -2.+2.*x[3], Jac[3, 4] = 0., Jac[4, 1] = 0., Jac[4, 2] = 0., Jac[4, 3] = 0., Jac[4, 4] = -2.+2.*x[4]]

 

irformJ := StatSeq(seq(Assign(op(i)),i=eqsJ)):

prcconsJ:= codegen[intrep2maple](Proc(Parameters(x::Array,Jac::Array),irformJ)):

J3:=Compiler:-Compile(prcconsJ):


 Following linear sovler algorithm is basiclly from dsolve/numeric/SC

.

s3:=proc(n::posint, A::Array( datatype = float[ 8 ] ) , ip::Array( datatype = integer[ 4 ] ),b::Array( datatype = float[ 8 ] ) )
local i::integer, j::integer, k::integer, m::integer, t::float;
      ip[n] := 1;
      for k to n-1 do
        m := k;
        for i from k+1 to n do
          if abs(A[m,k]) < abs(A[i,k]) then
            m := i
           end if
         end do;
        ip[k] := m;
        if m <> k then
          ip[n] := -ip[n]
         end if;
       t := A[m,k];
       A[m,k] := A[k,k];
       A[k,k] := t;
       if t = 0 then
         ip[n] := 0;
         #return ip
         end if;
       for i from k+1 to n do
         A[i,k] := -A[i,k]/t
         end do;
       for j from k+1 to n do
         t := A[m,j];
         A[m,j] := A[k,j];
         A[k,j] := t;
         if t <> 0 then
           for i from k+1 to n do
             A[i,j] := A[i,j]+A[i,k]*t
             end do
           end if
         end do
       end do;
     if A[n,n] = 0 then
       ip[n] := 0
       end if;
     #ip[n]
      if ip[n] = 0 then
        error "The matrix A is numerically singular"
       end if;
      if 1 < n then
        for k to n-1 do
          m := ip[k];
          t := b[m];
          b[m] := b[k];
          b[k] := t;
          for i from k+1 to n do
           b[i] := b[i]+A[i,k]*t
           end do
         end do;
       for k from n by -1 to 2 do
         b[k] := b[k]/A[k,k];
         t := -b[k];
        for i to k-1 do
           b[i] := b[i]+A[i,k]*t
           end do
         end do
       end if;
     b[1] := b[1]/A[1,1];
return 0.0;
end proc:

s3c:=Compiler:-Compile(s3):

Newton:=proc(x::Array(datatype=float[8]),
             f0::Array(datatype=float[8]),
             j0::Array(datatype=float[8]),
             N::integer,
             ip::Array(datatype=integer[4]))
global f3,J3,s3c;
local err::float[8],i::integer,k::integer,c::integer;
err:=10.0;
c:=0;
while err>1e-8 and c<100 do
  f3(x,f0);
  J3(x,j0);
  for i from 1 to N do
    for k from 1 to N do
      j0[i,k]:=-1.0*j0[i,k];
    end do;
  end do;
  s3c(N,j0,ip,f0);
  for i from 1 to N do
    x[i] := x[i]+f0[i];
  end do;
  f0[1] := abs(f0[1]);
  err := abs(f0[1]);
  for i from 2 to N do
    f0[i] := abs(f0[i]);
    if f0[i]>err then
      err := f0[i];
    end if;
  end do;
  err:=err/N;
  c:=c+1;
end do:
err;
end proc:

Newtonc := Compiler:-Compile(Newton):

Warning, the function names {J3, f3, s3c} are not recognized in the target language

x0:=Array(1..N,datatype=float[8]):
f0:=Array(1..N,datatype=float[8]):
j0:=Array(1..N,1..N,datatype=float[8]):
ip:=Array(1..N,datatype=integer[4]):

ArrayTools:-Fill(0.5,x0):
ArrayTools:-Fill(0.0,f0):
ArrayTools:-Fill(0.0,j0):
ArrayTools:-Fill(1,ip):
x0; f0; j0, ip;

t11:=time():
for ii from 1 to 10000 do
  Newtonc(x0,f0,j0,N,ip); # compiled
end do:
time()-t11;

x0, f0, j0, ip;

Array([.900000000000000, .132532421355126, 0.566018867943396e-1, 0.314314686094742e-1])

Array([0., 0., 0., 0.358203559119632e-17])

Array(1..4, 1..4, {(1, 1) = .20000000000000018, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (2, 1) = .0, (2, 2) = 1.7349351572897471, (2, 3) = -.0, (2, 4) = -.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 1.8867962264113207, (3, 4) = -.0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = 1.9371370627810516}, datatype = float[8]), Array(1..4, {(1) = 1, (2) = 2, (3) = 3, (4) = 1}, datatype = integer[4])

.430

Array(1..4, {(1) = .8999999999999999, (2) = .13253242135512638, (3) = 0.5660188679433962e-1, (4) = 0.31431468609474184e-1}, datatype = float[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = 0.3582035591196321e-17}, datatype = float[8]), Array(1..4, 1..4, {(1, 1) = .20000000000000018, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (2, 1) = .0, (2, 2) = 1.7349351572897471, (2, 3) = -.0, (2, 4) = -.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 1.8867962264113207, (3, 4) = -.0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = 1.9371370627810516}, datatype = float[8]), Array(1..4, {(1) = 1, (2) = 2, (3) = 3, (4) = 1}, datatype = integer[4])

ArrayTools:-Fill(0.5,x0):
ArrayTools:-Fill(0.0,f0):
ArrayTools:-Fill(0.0,j0):
ArrayTools:-Fill(1,ip):
x0; f0; j0, ip;

t11:=time():
for ii from 1 to 10000 do
  Newton(x0,f0,j0,N,ip); #uncompiled
end do:
time()-t11;

x0, f0, j0, ip;

Array([.900000000000000, .132532421355126, 0.566018867943396e-1, 0.314314686094742e-1])

Array([0., 0., 0., 0.358203559119632e-17])

Array(1..4, 1..4, {(1, 1) = .20000000000000018, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (2, 1) = .0, (2, 2) = 1.7349351572897471, (2, 3) = -.0, (2, 4) = -.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 1.8867962264113207, (3, 4) = -.0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = 1.9371370627810516}, datatype = float[8]), Array(1..4, {(1) = 1, (2) = 2, (3) = 3, (4) = 1}, datatype = integer[4])

1.380

Array(1..4, {(1) = .8999999999999999, (2) = .13253242135512638, (3) = 0.5660188679433962e-1, (4) = 0.31431468609474184e-1}, datatype = float[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = 0.3582035591196321e-17}, datatype = float[8]), Array(1..4, 1..4, {(1, 1) = .20000000000000018, (1, 2) = -.0, (1, 3) = -.0, (1, 4) = -.0, (2, 1) = .0, (2, 2) = 1.7349351572897471, (2, 3) = -.0, (2, 4) = -.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 1.8867962264113207, (3, 4) = -.0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = 1.9371370627810516}, datatype = float[8]), Array(1..4, {(1) = 1, (2) = 2, (3) = 3, (4) = 1}, datatype = integer[4])

eval(f,x=x0);

Vector[row](4, {(1) = 0., (2) = 0., (3) = 0., (4) = -0.6938893904e-17})

 


Download NewtonRaphsoncompile_1.mw

acer

something

May 05 2015 acer 11072

I made some changes. Mostly I had to guess at things.

For example I guessed that the loop to compute tx in the last part was intended as some sort of bisection thingy.

You don't need to load evalf using with. And it makes no sense to try that. Other packages like plots are loaded using with, not whit as you misspelled it in some places.

The lower end of a range passed to fsolve, in the first code part, needed to be less than 4.

It's slightly crazy that your code is loading the deprecated linalg package. I guess you saw that in a book or something, which is a shame.

The procedure dsnumsort is very poor. Did that come from a book too? It's methodology is awful. It's not surprising that you got confused trying to utilize it. I guessed at your purpose and utilized 2-argument eval instead. (How about using output=listprocedure for dsolve numeric, and perhaps make life easier?)

M11MAS94modif.mw

acer

another way

May 04 2015 acer 11072

Here's another way (where evalf/RoofOf will use fsolve, under the hood...), which happens to work for this particular example.

c := x+sin(x)+ln(t):
F := x+x^2:
plot(subs(x=solve(c,x),F),t=0..10,view=-1..6);

I have deliberately gone for simplicity here, rather than efficiency. Other ways are faster. But this is simple.

acer

indexing

May 01 2015 acer 11072

Yes, you can replace those two calls to ArrayTools:-Copy with a single indexed assignment.

restart;

A := LinearAlgebra:-RandomVector[row](10):
B := Vector[row](10):

B[[2,3,8,9]] := A[[2,3,8,9]]:

B;

                     [0, 27, 8, 0, 0, 0, 0, 92, -31, 0]

acer

something

April 30 2015 acer 11072
restart:

s := { x>4, y>4 }:

sbar := solve( map( `not`, s ) );

                          sbar := {x <= 4, y <= 4}

convert(s, RealRange);

       {x::(RealRange(Open(4), infinity)), y::(RealRange(Open(4), infinity))}

convert( sbar, RealRange );

             {x::(RealRange(-infinity, 4)), y::(RealRange(-infinity, 4))}

list

April 29 2015 acer 11072

Is is a small list of precomputed primes, used for quick look-up.

restart:
kernelopts(opaquemodules=false):
isprime:-special_primes;

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 

  73, 79, 83, 89, 97, 979073763359, 987715808641, 1203502041253, 987709929361, 

  987714327601, 902164268009, 18945002166578281, 987723662641, 987674160001]

acer

ssystem

April 26 2015 acer 11072

The ssystem command can execute a command on the host OS, and is available across platforms and Maple interfaces.

acer

isolate

April 15 2015 acer 11072

How about using `isolate` instead of `solve`?

rhs(isolate(Eq1, diff(y1(t), t, t)));

acer

prettyprint?

April 14 2015 acer 11072

So you working in Matlab, with Maple as the engine for its Symbolic Toolbox, is that right?

And you want results from Maple (invoked from within Matlab) to be printed in the Matlab interface as something which then cut&pastes as 1D Maple Notation input? Is that right?

Perhaps it might help to first issue the (Maple command) interface(prettyprint=0)

acer

bug

April 14 2015 acer 11072

Looks like a bug.

Try it instead as,

LinearAlgebra:-Eigenvectors( convert(A,rational) );

I will submit a bug report.

acer

goals

April 10 2015 acer 11072

Please describe in as much detail as possible how you want color to be defined.

The code you're shown has procedure mandelbrot which currently returns a single value given input point (x,y). And plot3d is using that to select a hue (as in H=hue from the HSV color space).

With just a single returned value m (per x,y input) it's natural to either produce a grayscale shading value or a hue value. Or you could use m for the intensity (the V in HSV) as well as for the hue.

Another value that you *could* make use of is the final z, for each (x,y) input. Since z is complex its real absolute value abs(z) is another piece is information that can be combined with the final m value to produce a color. With two values (m and abs(z) say) you could produce some attractive color shadings by generating three numbers for the RGB color space. Or you could produce values in a 3-element color space by doing modular arithmetic on either a single or a double-valued return.

Note that I have made no special effort here for efficiency. In particular the 2-value-return version mandelbrot2 below is not evalhf'able and so is slow. It could be made faster.

I didn't understand why you'd want a style=point, so I changed it. See also the last example on the densityplot command's help page. Other related links: herehere, here, and (distantly, showing how images and plot GRID are interchangeable) here and here.

The plots below look much brighter and less grainy in Maple itself.

restart:

mandelbrot := proc(x, y)

     local c, z, m;

     c := evalf(x+y*I);

     z := c;

     for m from 0 to 30 while abs(z) < 2.0 do

        z := z^2+c

        od;

        m;

     end:

plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[200, 200],
       scaling=constrained, style=patchnogrid, lightmodel=none,
       color=[((x,y)->mandelbrot(x,y)/30.0),
              1,
              1,
              colortype=HSV]);

plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[200, 200],
       scaling=constrained, style=patchnogrid, lightmodel=none,
       color=[((x,y)->mandelbrot(x,y)/30.0),
              1,
              ((x,y)->30-mandelbrot(x,y)),
              colortype=HSV]);

p1,p2,p3 := 23,5,2:

plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[200, 200],
       scaling=constrained, style=patchnogrid, lightmodel=none,
       color=[((x,y)->irem(mandelbrot(x,y),p1)),
              ((x,y)->irem(mandelbrot(x,y),p2)),
              ((x,y)->irem(mandelbrot(x,y),p3)),
              colortype=RGB]);

mandelbrot2 := proc(x, y)
     option remember,system, hfloat;

     local c, z, m;

     c := x+y*I;

     z := c;

     for m from 0 to 30 while abs(z) < 2.0 do

        z := z^2+c

        od;

        m, abs(z); # no longer evalhfable, as is

     end:

plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[200, 200],
       scaling=constrained, style=patchnogrid, lightmodel=none,
       color=[((x,y)->mandelbrot2(x,y)[1]/30.0),
              1,
              ((x,y)->mandelbrot2(x,y)[2]),
              colortype=HSV]);

 


Download naive_mandelbrot.mw

acer

like this?

April 09 2015 acer 11072

I am having trouble understanding what you want, sorry. Perhaps something like this?

(It's a bug in this server that gridlines are shown in the first plot below.)

restart;

R:=2:

plots:-display(
   plots:-polarplot(R,theta=0..2*Pi,color=red),
   plots:-polarplot(R+cos(theta),theta=0..2*Pi,color=blue),
   gridlines=false);

plots:-display(
   plot(R,theta=0..2*Pi,color=red,coords=polar),
   plot(R+cos(theta),theta=0..2*Pi,color=blue,coords=polar),
   gridlines=false,scaling=constrained);

 


Download costhetafun.mw

acer

1 2 3 4 5 6 7 Last Page 1 of 109