acer

32587 Reputation

29 Badges

20 years, 37 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer


restart:

with(plots):

nx:=20: tmax:=50: T1:=1: T2:=10: L:=1: k:=1: rho:=1: cp:=1:
chi:=k/rho/cp: h:=L/(nx-1): t:=1e-3:

T:=Array(0..nx,0..tmax,datatype=float[8]):
for k from 0 to nx do T[k,0]:=T1 od:

for w from 1 to tmax do

  T[0,w]:=T1: T[nx,w]:=T2:

  for q from 1 to nx-1 do

    T[q,w]:=T[q,w-1]+chi*t/h^2*(T[q+1,w-1]+T[q-1,w-1]-2*T[q,w-1]);

  od:

od:

surfdata(Matrix(T), 0..nx, 0..tmax, style=surface,

         dimension=2, labels=["x", "t"],
         #dimension=3, labels=["x", "t", "T"],

         gridsize=[100,100], # interpolate, for smoothness

         colorscheme=["zgradient",
                      ["Indigo","Blue","Cyan","Green","Yellow","Orange","Red"]]);


Download tempdens.mw

acer

See the help topic updates/Maple2016/compatibility in your Maple 2016.


restart;

kernelopts(version);

`Maple 2016.0, X86 64 WINDOWS, Feb 16 2016, Build ID 1113130`

S:=convert(1/(3*x+2),FormalPowerSeries);

Sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)

_EnvFormal;

_EnvFormal

value(S);

sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)

sum(op(S), formal);

1/(3*x+2)

restart;

S:=convert(1/(3*x+2),FormalPowerSeries):

_EnvFormal:=true:

value(S);

1/(3*x+2)

restart;

S:=convert(1/(3*x+2),FormalPowerSeries):

_EnvFormal:=false:

value(S);

sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)

value(S) assuming x<2/3, x>-2/3;

1/(3*x+2)

 


Download formalsum2016.mw

acer

(It turns out that I mostly just had to look for the file where I'd done pretty much this, earlier.)

If the IterativeMaps commands also give you grief you could also implement it (slower) using plot3d or densityplot.

As mentioned in the code below, increasing `maxiter` will reduce the areas of nonconvergence (black areas) for such univariate polynomials.

I used Maple 2015.2 on 64bit Linux.

restart:

with(IterativeMaps):
with(ImageTools):

ee := z^5 - 1;

z^5-1

## For ee a polynomial we can compute all the roots.
## This is used to specify the colors, below.
rts:=[fsolve(ee,z,complex)]:
evalf[4](%);

[-.8090-.5878*I, -.8090+.5878*I, .3090-.9511*I, .3090+.9511*I, 1.]

minr,maxr := (min,max)(Re~(rts)):
mini,maxi := (min,max)(Im~(rts)):

dee := diff(ee, z):
Z := eval( z - ee/dee, z=zr+zi*I):

fzi := evalc( Im( Z ) ):

fzr := evalc( Re( subs(zi=zit, Z) ) ):

N := 400;
maxiter := 40; # increasing this will reduce nonconvergence
h,w := N, floor(N*(maxr-minr)/(maxi-mini));

400

40

400, 380

newt := Escape( height=h, width=w,
                       [zi, zr, zit],
                       [fzi, fzr, zi],
                       [y, x, y], evalc(abs(eval(ee,[z=zr+zit*I])))^2<1e-10,
                       2*minr, 2*maxr, 2*mini, 2*maxi,
                       iterations = maxiter,
                       redexpression = zr, greenexpression = zit, blueexpression=i ):

## View the iteration count layer
#Embed([ FitIntensity(newt[..,..,3]) ]);

## By default use full saturation and full intensity
img:=Array(1..h,1..w,1..3,(i,j,k)->1.0,datatype=float[8]):

## Optional, set intensity to zero for non-convergence
img[..,..,3]:=Threshold(newt[..,..,3],1,low=maxiter+1):
img[..,..,3]:=Threshold(img[..,..,3],maxiter-1,high=0,low=1):

## Optional, scale the saturation by the time to converge
img[..,..,2]:=1.0-~FitIntensity(newt[..,..,3]):

## Done inefficiently, set the hue by position in the list of
## all roots
G:=Array(zip((a,b)->min[index](map[evalhf](abs,rts-~(a+b*I))),
             newt[..,..,1],newt[..,..,2]),datatype=float[8],order=C_order):
img[..,..,1]:=240*FitIntensity(G):

## Convert from HSV to something viewable
final := HSVtoRGB(img):

Embed(final);

NULL

Download newtbasin0.mw

 

And here it is for N=600 and maxiter=200. Using a logarithm for adjusting the saturation layer might be a way to get a nicer impression of the convergence rate.

 

acer

c:=INTERVAL(1.41421356167 .. 1.41421356308):

uo,low:=op(op(c));                          

                           uo, low := 1.41421356167, 1.41421356308

Or,

up,low:=op([1,..],c);

                           up, low := 1.41421356167, 1.41421356308

acer

You can use the add and the if command together.

Note the single backticks (left-quotes) around the word if.

V := Vector[row]( 10, i -> i^2 );

                 V := [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]

add( `if`( v>30, v, 0), v=V );

                                     330

add( map( v -> `if`( v>30, v, 0 ), V ) );

                                     330

add( piecewise( v>30, v, 0 ), v=V );

                                     330

M := Matrix( 5, 5, (i,j) -> (i+j) );

                                 [2  3  4  5   6]
                                 [              ]
                                 [3  4  5  6   7]
                                 [              ]
                            M := [4  5  6  7   8]
                                 [              ]
                                 [5  6  7  8   9]
                                 [              ]
                                 [6  7  8  9  10]

add( `if`( M[i,3] > 5, M[i,3], 0), i=1..5 );

                                     21

add( `if`( v > 5, v, 0), v = M[..,3] );

                                     21

Is that the kind of thing you mean?

acer

 

restart;

#
# One way to consider the problem is to notice that the subterm
# 18*s^2+32 happens to equal  9*s^2  +  9*s^2+32  where those
# two parts are the M^2 and N^2 of the desired (M+N)^2 target.
#
# We can forcibly separate those enough to allow factor or CompleteSquare
# to obtain the desired form. But how to automate the split?
#

a := 6*s*sqrt(9*s^2 + 32) + 18*s^2 + 32;

6*s*(9*s^2+32)^(1/2)+18*s^2+32

t := indets(a, `^`(anything, 1/2))[1];

(9*s^2+32)^(1/2)

ft := freeze(t);

`freeze/R0`

subs(t = ft, a - t^2 + ft^2); # a rather manual step

9*s^2+6*s*`freeze/R0`+`freeze/R0`^2

thaw(factor(%));

(3*s+(9*s^2+32)^(1/2))^2

# But how do we know how to divvy it up? How to find A=1?
subs(t = ft, a + A*(-t^2 + ft^2));
Student:-Precalculus:-CompleteSquare(%, ft);
subs(S=s^2, collect(algsubs(s^2=S, %), [S], u->simplify(u, size)));
thaw( eval(%, A=1) );

6*s*`freeze/R0`+18*s^2+32+A*(-9*s^2+`freeze/R0`^2-32)

A*(`freeze/R0`+3*s/A)^2+18*s^2+32+A*(-9*s^2-32)-9*s^2/A

-9*(A-1)^2*s^2/A+A*(`freeze/R0`+3*s/A)^2-32*A+32

(3*s+(9*s^2+32)^(1/2))^2

restart;
# Once again, if you know in advance how to split terms, it's easy.
# How to automate it?
eq1 := 6*s*sqrt(9*s^2 + 32)=2*M*N:
eq2 := 18*s^2 + 32=M^2+N^2:
map(factor, eq1+eq2);
eval(%, solve( {eq1,eq2}, {M,N}, explicit )[1]);

6*s*(9*s^2+32)^(1/2)+18*s^2+32 = (M+N)^2

6*s*(9*s^2+32)^(1/2)+18*s^2+32 = ((9*s^2+32)^(1/2)+3*s)^2

 

 

Download factorfun.mw

acer

If the following code is saved to a Maple initialization file then a new entry for calling the RationalFunctionPlot command should appear in the right context-menu, under the existing "Plots" submenu, if the expression is of type `ratpoly`.

You can of course adjust this. The `caption` might not be set to the empty string. Or several such new menu items could be added, in their own subsubmenu, to provided various common choices for the command's options.

:-__AugmentCM := :-ContextMenu:-CurrentContext:-Copy():
:-__AugmentCM:-Entries:-Add("Rational Function",
            ":-Student:-Precalculus:-RationalFunctionPlot(%EXPR,':-caption'=\"\")",
            ':-ratpoly',
            ':-operator'=:-Typesetting:-mover(Typesetting:-mo("→"),
                         :-Typesetting:-mi("Rational Function Plot")),
            ':-category'="Category 1",
            ':-helpstring'="Rational Function Plot",
            ':-submenu'=["Plots"]
           ):
:-ContextMenu:-Install(:-__AugmentCM);

Or you could adjust this code to invoke instead the Student:-Precalculus:-RationalFunctionTutor command. This command is already available via the main menubar choice Tools->Tutors->Precalculus but invoking it from there won't automatically pass in the desired expression. But if this command is added to the context-menu then the right-clicked expression will be passed in and used.

See the help page with topic worksheet,reference,initialization for documentation on how to create an initialization file.

If you want all the students to get this enhancement automatically for a network installation of Maple then you just add it to the (possibly new) file <MAPLE>/lib/maple.ini where <MAPLE> is the location of the Maple installation. Or do it for each machine, if Maple is installed locally on machines. Test it first, of course.

acer

h := x -> f(g(x)):

D(h)(0);

                             D(f)(g(0)) D(g)(0)

eval( D(h)(0), D(g)(0)=0 );

                                      0

Even if you prefer `diff` over `D` it still evaluates to zero, for the same reason.

eval( eval(diff(h(x),x),x=0), eval(diff(g(x),x),x=0)=0 );

                                      0

acer

 

restart;

f := (125*x^3+525*x^2*y+735*x*y^2+343*y^3-36*z^2)
     /(-216*z^3+25*x^2+70*x*y+49*y^2);

(125*x^3+525*x^2*y+735*x*y^2+343*y^3-36*z^2)/(-216*z^3+25*x^2+70*x*y+49*y^2)

collect(f,z,factor);

(-36*z^2+(5*x+7*y)^3)/(-216*z^3+(5*x+7*y)^2)

Download collectfactor.mw

acer

This is covered in the subsubsubsection "2-dimensional case"->"Interpolation"->"Uniform data grid"->"Interpolating procedure" of the help page available via ?examples,Interpolation_and_Smoothing .

The basic methodology is simple, see `f` below. For that I'd suggest using Vectors/Array/Matrices and square-bracket indexing and not lists and `op`.

Using float[8] datatype helps with resources. There is overhead in getting a pair of scalars into the indexable containers (Arrays, lists, etc) that ArrayInterpolation accepts. It's designed to be fast for dealing with large amounts of interpolating points all at once. Repeated use to process individual points will be less efficient than processing them all at once. But perhaps you won't know up front what they all are. There is an efficiency section in that example worksheet for the 1D interpolation case's efficiency, and I follow that is the more efficienct `F` below.


restart;

x:=<0,1,2>:
y:=<0,1,2>:
data:=<<0,2,4>|<1,4,16>|<1,8,64>>:

f:=(a,b)->CurveFitting:-ArrayInterpolation([x,y], data,
              Array(1..1, 1..1, 1..2, [[[a,b]]]), method=spline)[1,1]:

f(0.5,0.5);

HFloat(1.5537109375)

#plot3d(f, 0..2, 0..2, style=surface);

#dataplot(x, y, data, style=surface);

#plots:-surfdata(Array(data,datatype=float[8]), 0..2, 0..2, source=regular, style=surface);

#plots:-matrixplot(data, labels=["","",""], style=surface,
#                  axis[1]=[tickmarks=[]], axis[2]=[tickmarks=[]]);

 

#
# The following `F` is more efficient than `f` above.
#
F := module()
export ModuleApply;
local inArray, outArray, hwxdata, hwydata, hwdata;
   inArray:=Array(1..1,1..1,1..2,':-datatype'=':-float[8]');
   outArray:=Array(1..1,1..1,':-datatype'=':-float[8]');
   hwxdata,hwydata := Vector(x,':-datatype'=':-float[8]'),
                      Vector(y,':-datatype'=':-float[8]'):
   hwdata:=Matrix(data,':-datatype'=':-float[8]'):
   ModuleApply:=proc(a, b,
                    {method::identical(spline,linear,cubic):=':-spline'})
     if not ( a::numeric and b::numeric ) then
       return 'F'(args);
     end if;
     UseHardwareFloats:=true;
     inArray[1,1,1]:=a;
     inArray[1,1,2]:=b;
     CurveFitting:-ArrayInterpolation([hwxdata,hwydata], hwdata,
                                      inArray, ':-container'=outArray,
                                      ':-method'=method, _rest);
      outArray[1,1];
   end proc;
end module:

f(0.5,0.5), F(0.5,0.5);
f(0.3, 0.45), F(0.3, 0.45);
F(s,t);

HFloat(1.5537109375), HFloat(1.5537109375)

HFloat(1.1073818203124999), HFloat(1.1073818203124999)

F(s, t)

CodeTools:-Usage( plot3d(f, 0..2, 0..2) ):

memory used=78.18MiB, alloc change=32.00MiB, cpu time=4.34s, real time=4.34s, gc time=93.60ms

CodeTools:-Usage( plot3d(F, 0..2, 0..2) ):

memory used=33.87MiB, alloc change=0 bytes, cpu time=1.82s, real time=1.83s, gc time=46.80ms

#plot3d(F, 0..2, 0..2, style=surface);

#plot3d(F(s,t,method=linear), s=0..2, t=0..2, style=surface);

 


Download 2dinterp.mw

acer

Apart from increasing the grid sizing to get a smoother 3D plot (as Doug suggested), you can also make the 2D contourplot look a little better by using a different choice of custom contour values.

Download sz2012_4m_K.mw

[edited] The 2D colorbar facsimile (faked, as a legend) can look a little nicer if the `.` used to force the nonempty strings (color blobs) are also colored.

sz2012_4m_KK.mw

The GUI seems unwilling to typeset colored 2D Math if the mathcolor (foreground color) and the mathbackground (background color) are exactly the same, on my Linux at least. A small adjustment of either, so they are not identical, seems to work. Proper care should be taken to make the adjusted RGB24 color value be valid with entries all in 0-255 (which I didn't do). A further possible refinement is that intsead of atomic-identifiers like `#mrow(...)` one could use Typsettting:-mrow, etc, so that size=posint qualifier could be used to also modify the typeset size. (I didn't do that either.)

 

restart;

PP:=sin(x)*y^2:

(c1,c2):=ColorTools:-Color("#7700ff"),ColorTools:-Color("#33ff66"):

contlist:=[seq(i, i=-1..1, .125)]:

numcont:=nops(contlist):

clist:=[seq(ColorTools:-Color([c1[]]+(i-1)*([c2[]]-[c1[]])/(numcont-1)), i=1..numcont)]:

P1:=plots:-contourplot(PP, x=-Pi..Pi, y=-1..0, filledregions, grid=[100,100],
                     coloring=[clist[1], clist[-1]], contours=contlist):

P2:=plots:-display(
    seq(plot([[1,0]], color=white, axes=none,
           legend=typeset(nprintf(cat("#mrow(mi(\".      .\",mathbackground=%a,mathcolor=%a)",
                                      `if`(irem(i,2)=1,
                                           sprintf(",mi(\"gt;\"),mn(\"% .2f\")",
                                                   contlist[i]),NULL),")"),
                                  ColorTools:-RGB24ToHex(ColorTools:-RGBToRGB24([clist[i+1][]])),
                                  ColorTools:-RGB24ToHex(ColorTools:-RGBToRGB24([clist[i+1][]])
                                                         -[0,0,1])))),
      i=nops(clist)-1 .. 1, -1),legendstyle=[location=right]):

plots:-display(P1, P2, legendstyle=[location=right], size=[680,500], axes=box,
              axis[1]=[tickmarks=piticks], labels=[``,``]);

 

Download altcolorbar.mw

This legend-as-colorbar trick could be made into a convenient procedure. It doesn't have to be used with just contourplots. It might also be used with densityplots, or 2D pointplots, or contour plots implemented via implicitplot calls.

 

acer

The problem is that when you save a Units:-Unit beast to a .mla archive and then access it in a new session the attributes (that commands in Units:-Standard expect to see on it) have been stripped off. That makes the command error out.

This can be illustrated by example:

restart:

A:=Units:-Unit(m):
lprint(A);

Units:-Unit(m)

attributes(A);

Units:-UnitStruct(1, metre, contexts:-SI), inert

Units:-Standard:-`*`(A, 2.5, A);

2.5*Units:-Unit(('m')^2)

fn:=cat(kernelopts(homedir),"/My Documents/foobar.mla"):

try
  LibraryTools:-Create(fn);
catch "there is already an archive":
end try;

LibraryTools:-Save('A', fn);

Download unitmla1.mw

 

Now let's restart and access that saved value. The same call to Units:-Standard:-`*` now emits an error message. This happens because the internal code of that command is expecting those attributes, as printed above.

 

restart;

libname:=cat(kernelopts(homedir),"/My Documents/foobar.mla"), libname:

lprint(A);

Units:-Unit(m)

attributes(A); # NULL

Units:-Standard:-`*`(A, 2.5, A); # error due to lack of expected attributes

Error, (in Units:-Standard:-*) invalid subscript selector

Download unitmla2.mw

 

I'm not sure what the most practical workaround to this problem might be. The following reinstantiates the Units:-Unit call and produces the value anew, with the attributes. I don't know whether other, related problems might arise though.

 

restart;

libname:=cat(kernelopts(homedir),"/My Documents/foobar.mla"), libname:

lprint(A);

Units:-Unit(m)

fixit:=(u)->subsindets(u,':-specfunc(anything,Units:-Unit)',convert,`global`):

A:=fixit(A):

attributes(A);

Units:-UnitStruct(1, metre, contexts:-SI), inert

Units:-Standard:-`*`(A, 2.5, A);

2.5*Units:-Unit(('m')^2)

Download unitmla3.mw

 

This behaviour appears for me in Maple 13, 17, and 2015 (which are all I tested). I shall submit a bug report.


acer


fo:=f(1/x);

f(1/x)

diff(fo,x); # the symbolic derivative

-(D(f))(1/x)/x^2

convert(%,diff); # an alternate representation

-(eval(diff(f(t1), t1), {t1 = 1/x}))/x^2


Download diffat.mw

acer

There is a command Statitistics:-ExponentialFit which allows you to weight the data (as applied to a transformed model -- see the help page).

You may note that I am not supplying initial points here.

You should probably look at some of the other pieces of information that the `output` option of this command can supply to you. A central question which you ought to be asking is: how are you going to judge one fit as being better than another? Only you can decide that, and until you do it isn't very sensible for the rest of us to use any metric to make such judgement.

restart;

year := [1850,1860,1870,1880,1890,1900,1910,1920,1930,1940,
         1950,1960,1970,1980,1990,2000,2010]:

population :=[4668,9070,17375,27985,37249,63786,115693,
              186667,359328,528961,806701,1243158,1741912,
              2049527,2818199,3400578,4092459]:

sol:=Statistics:-ExponentialFit(year, population, x);

                                     -31                          
           sol := 1.08291520188752 10    exp(0.0433611679840687 x)

plots:-display(plot(sol,x=year[1]..year[-1]),
               plot(<<year>|<population>>,style=point));

W := map(y->abs(1850-y)^3,year):

solW:=Statistics:-ExponentialFit(year, population, x, weights=W);

                                     -22                          
          solW := 6.74209319468292 10    exp(0.0319134519471405 x)

plots:-display(plot(solW,x=year[1]..year[-1]),
               plot(<<year>|<population>>,style=point));

 

Download expfit.mw

acer

Using the same code as I did in your duplicate post:


restart;

eq := (-sin(p(t)+g(t))*cos(a(t))-sin(b(t))*sin(a(t))*cos(p(t)+g(t)))*l2[1]+((-sin(p(t)+g(t))*cos(a(t))-sin(b(t))*sin(a(t))*cos(p(t)+g(t)))*cos(th(t))+(-cos(p(t)+g(t))*cos(a(t))+sin(a(t))*sin(b(t))*sin(p(t)+g(t)))*sin(th(t)))*l3[1]+(-sin(a(t))*sin(g(t))*sin(b(t))+cos(a(t))*cos(g(t)))*xb[1]+sin(a(t))*cos(b(t))*yb[1]+(sin(a(t))*sin(b(t))*cos(g(t))+cos(a(t))*sin(g(t)))*zb[1]+x(t)-xp(t) = 0;

(-sin(p(t)+g(t))*cos(a(t))-sin(b(t))*sin(a(t))*cos(p(t)+g(t)))*l2[1]+((-sin(p(t)+g(t))*cos(a(t))-sin(b(t))*sin(a(t))*cos(p(t)+g(t)))*cos(th(t))+(-cos(p(t)+g(t))*cos(a(t))+sin(a(t))*sin(b(t))*sin(p(t)+g(t)))*sin(th(t)))*l3[1]+(-sin(a(t))*sin(g(t))*sin(b(t))+cos(a(t))*cos(g(t)))*xb[1]+sin(a(t))*cos(b(t))*yb[1]+(sin(a(t))*sin(b(t))*cos(g(t))+cos(a(t))*sin(g(t)))*zb[1]+x(t)-xp(t) = 0

length(eq), MmaTranslator:-Mma:-LeafCount(eq), `simplify/size/size`(eq);

780, 148, 321

ans1 := collect(expand(eq), [l2[1],l3[1]],
                u->simplify(frontend(combine, [u],
                                     [{`+`,`*`},
                                      {seq('cos(x(t)),sin(x(t))',x=[p,g,th])}]),
                            size));

(-sin(p(t)+g(t))*cos(a(t))-sin(b(t))*sin(a(t))*cos(p(t)+g(t)))*l2[1]+(-cos(a(t))*sin(g(t)+th(t)+p(t))-sin(b(t))*sin(a(t))*cos(g(t)+th(t)+p(t)))*l3[1]+(cos(g(t))*sin(b(t))*zb[1]-sin(b(t))*sin(g(t))*xb[1]+yb[1]*cos(b(t)))*sin(a(t))+(cos(g(t))*xb[1]+sin(g(t))*zb[1])*cos(a(t))+x(t)-xp(t) = 0

length(ans1), MmaTranslator:-Mma:-LeafCount(ans1), `simplify/size/size`(ans1);

590, 111, 250

ans2 := simplify(frontend(combine, [expand(eq)],
                          [{`+`,`*`},
                           {seq('cos(x(t)),sin(x(t))',x=[p,g,th])}]),
                 size);

-l3[1]*sin(b(t))*sin(a(t))*cos(g(t)+th(t)+p(t))-l3[1]*cos(a(t))*sin(g(t)+th(t)+p(t))-l2[1]*sin(b(t))*sin(a(t))*cos(p(t)+g(t))-l2[1]*cos(a(t))*sin(p(t)+g(t))+((zb[1]*cos(g(t))-sin(g(t))*xb[1])*sin(b(t))+yb[1]*cos(b(t)))*sin(a(t))+(cos(g(t))*xb[1]+sin(g(t))*zb[1])*cos(a(t))-xp(t)+x(t) = 0

length(ans2), MmaTranslator:-Mma:-LeafCount(ans2), `simplify/size/size`(ans2);

586, 108, 241

simplify(eq - ans1), simplify(eq - ans2);

0 = 0, 0 = 0

 


Download tringsimpv.mw

acer

First 213 214 215 216 217 218 219 Last Page 215 of 338