acer

32348 Reputation

29 Badges

19 years, 329 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

If you have assigned your original inert integral of a piecewise to the name K, then you could do,

value(K) assuming continuous:
simplify(%) assuming T2>T1;

acer

Inspired by Rouben's very nice explanation, the result can be demonstrated via computation alone.

(The introduced values assigned to `a` and `bT` are inspired by RR's visual and verbal explanation, yet while not a recipe to solve all similar problems the following does comprise a fully computational result, FWIW.)


restart:

J := Int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1);

Int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1)

s := solve(op(1,J));

(1/2)*2^(4/5)

# J = J1 + c

J1,c := op(simplify(IntegrationTools:-Split(J,s))):
J = J1 + c;

Int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1) = Int(-x+(-x^5+1)^(1/5), x = 0 .. (1/2)*2^(4/5))+Int(x-(-x^5+1)^(1/5), x = (1/2)*2^(4/5) .. 1)

a := Int(s-x,x=0..s):

# J1 = a + b , hence J = a + b + c

b := IntegrationTools:-Combine(J1-a):
J1 = a + b;

Int(-x+(-x^5+1)^(1/5), x = 0 .. (1/2)*2^(4/5)) = Int((1/2)*2^(4/5)-x, x = 0 .. (1/2)*2^(4/5))+(1/2)*(Int(2*(-x^5+1)^(1/5)-2^(4/5), x = 0 .. (1/2)*2^(4/5)))

bT := Int((-x^5+1)^(1/5),x = 1/2*2^(4/5) .. 1);

# Demonstrate that bT = b , hence J = a + bT + c
# This shows that bT = b, though it might also serve
# to obtain bT from b, via a change of variables.

value(combine(bT-student[changevar](x=y^(1/5),

                   (student[changevar](y=(1-x^5),b,y)),x)));

Int((-x^5+1)^(1/5), x = (1/2)*2^(4/5) .. 1)

0

# Transform aT = a , hence J = aT + bT + c

aT := student[changevar](x=y,student[changevar](y=s-x,a,y),x);

Int(x, x = 0 .. (1/2)*2^(4/5))

J = aT + bT + c;
value(combine(%));

Int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1) = Int(x, x = 0 .. (1/2)*2^(4/5))+Int((-x^5+1)^(1/5), x = (1/2)*2^(4/5) .. 1)+Int(x-(-x^5+1)^(1/5), x = (1/2)*2^(4/5) .. 1)

int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1) = 1/2

 


Download rr.mw

acer

Using Maple 2015.1 and its Import command, letters can be imported from the STL 3D format.

restart:

#
# One (of many) STL format collections of the English alphabet:
#     http://www.stlfinder.com/model/all-alphabet-letters-a-z
#
basedir:=cat(kernelopts(homedir), "/My Documents/stl/All_Alphabet_Letters_A-Z"):

getletter:=L->Import(cat(basedir, "/Letter_", L, ".stl"),
                     orientation=[-90,0,0],title=""):

letters := [m,a,p,l,e]:

for letter in letters do
  if not assigned(AZ[letter]) then
    AZ[letter]:=getletter(letter);
  end if;
end do:

with(plottools): with(plots):

Rplot:=display(AZ[letters[1]],
               seq(translate(AZ[letters[i]], 40*(i-1)+15, 0, 0),
                   i=2..nops(letters))):

Rplot:=transform((x,y,z)->[z,x,y])(scale(Rplot, 5/6*2*Pi/200, 1/5, 1/20, [0,0,10])):

display(changecoords(Rplot,cylindrical), scaling=constrained, axes=none,
        orientation=[0,90,0], viewpoint=[circleright], frames=50);

acer

The first variant given on your cited link, with the writing on an opaque surface, is demonstrated below.

restart:
with(plots): with(ImageTools):
fn:=cat(kernelopts(homedir),"/mma.bmp"):
p:=display(textplot([0,1,"Mathematica",font=["times","roman",50]]),
           axes=none,scaling=constrained):
plotsetup(bmp,plotoutput=fn);
display(p);
plotsetup(default);
img:=Scale(Read(fn),1..250,1..250):
P3D:=plot3d(1, theta=0..2*Pi, h=0..Pi, image=img, coords=cylindrical,
            view=1.2..1.9, orientation=[0,90,0], axes=none,
            viewpoint=["circleright", frames=30]):
P3D;

Just like a usual animated Maple plot, you can single left-click on the plot, and then from the main menubar you could bump up the frame rate, and toggle loop, and play it.

If a non-rectilinear 3D projection effect is desired then a custom viewpoint formula or (more likely) a sequence of orientations is possible, instead of the simple circleright viewpoint.

acer

If the x and y (ie, the x_i) only appear in that way (as polynomial coefficients, say) then could you not simply evaluate the Matrix at [x_1=1,...x_n=1]?

M:=Matrix([[a*x, -b*y],[c*x^2, -d*y]]);

                                   [a x   -b y]
                                   [          ]
                              M := [   2      ]
                                   [c x   -d y]

eval(M, [x=1,y=1]);

                                   [a  -b]
                                   [     ]
                                   [c  -d]

acer

Use add rather than sum for such finite, definite summation. 

The problem you saw is that V[i] with an unassigned symbolic index i is not allowed for the kind of data structure associated with a DataTable.

For the sum command Maple evaluates the first argument up front, and baulks as V[i]. The add command has special evaluation rules (to delay the evaluation) and so Maple doesn't then evaluate V[i] until i gets a value.

Another way is to use unevaluation single right quotes, to delay the evaluation. Eg,

sum('V[i]', i=1..10);

acer

This works in at least 64bit Linux Maple 2015.0, 18.02, and 17.02,

ee := sin(-(1/6)*Pi+(1/2)*arccos(1/3)):

evala(expand(convert(ee,expln)));      

                                       1/2  1/2
                                      2    3
                                1/2 - ---------
                                          6

acer

Open a new blank Worksheet, and from the main menubar select Format -> Styles...

That should pop up a small window titled "Style Management". In its list box scroll down until you see "Maple Input". Select that. Then press click on the "Modify" button in that popup.

That should raise another popup window titled "Character Style" with an entry box filled in already with "Maple Input", and default "Courier New" grayed out. On the left of this second popup is a button "Color". Press it, and a color choose pops up. Select the bright red square. The press the "OK" button in that second popup, which dismisses it. Then press the "OK" button in the first popup, which dismisses it.

The 1D Maple Notation in execution groups and their prompts should now appear in bright red.

Now you need to make this blank worksheet your new default style. From the menubar select Format -> Manage Style Sets...

In the popup titled "Style Set Management" which appears click on the button "New Style Set". That raises another popup, titled "Choose Styles". In this popup clickon the "Unselect All" button to clear all the check boxes. Scroll down and toggle the check box "Maple Input". Then click on the "OK" button, which should make the file-manager popup appear, to allow you to save the style as a .mw file. Give it a name of your choice, and click on the "Save" button which dismisses the file-manager popup. You should now be back in the "Style Set Management" popup. Ckick the radio-button for "User-defined Style Set", and then click on the button "Browse" and find and choose the file which you previously saved. Then click on the "OK" button to finish.

At this point your new custom style (.mw file) sheet will get used for all new Worksheets, even if you close the whole GUI down and relaunch.

See here.

acer

In Maple 2015 you can programmatically insert a PlotComponent which contains the animation and then also set it playing programmatically.

When it's playing, you can still click on it with the mouse to get the usual menubar controls, so as to stop it, etc.

Here's a command for doing that. It needs at least Maple 2015.0 to execute the procedure.

Note that only one such insertion can be done per execution group.

Both the examples play the animation immediately upon insertion. The second one also loops (ie. replays continuously).

The MapleNet backend for rendering uploaded worksheets seems pretty old, and doesn't show the components. So below I've inserted then by hand into this posting. In the actual worksheet there is a border around each PlotComnent (and the code could be adjusted to repress that, too).

restart:

kernelopts(version);

`Maple 2015.1, X86 64 WINDOWS, May 14 2015, Build ID 1044754`

autoplay := proc( anim, {continuous::truefalse:=false} )
  local P, T;
  uses DocumentTools, DocumentTools:-Layout,
       DocumentTools:-Components;
  P := Plot(':-identity' = "Plot0", anim,
            ':-continuous'=continuous);
  T:=InsertContent(Worksheet(Group(Input(Textfield(P)))),
                   ':-output'=':-table');
  DocumentTools:-SetProperty(T["Plot0"],':-play',true);
  NULL:
end proc:

A := plots:-animate(plot,
                    [sin(a*x)/x^2,x=-Pi/2..Pi/2,
                     view=-30..30,gridlines=false],
                    a=1.0..50.0):

autoplay(A);

autoplay(A, continuous);

 

Download autoplay.mw

I have another, unfinished, fancier version which also has all the buttons to play/pause/etc right beside the PlotComponent. Not sure how much demand there'd be for that...

acer

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

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

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

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

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

First 225 226 227 228 229 230 231 Last Page 227 of 336