acer

32343 Reputation

29 Badges

19 years, 326 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Using Digits:=15 my 64bit Maple 2017.1 for Linux obtained a solution with a pretty good residual, in not so bad time.


 

restart;

kernelopts(version);

`Maple 2017.1, X86 64 LINUX, Jun 19 2017, Build ID 1238644`

Digits:=15:

# m:=2.5: N:=15: h:=0.29669:
m:=3: N:=17: h:=0.41600:

p := proc(x)
          c[-N-1]*x^2+1
     end proc:
dp := diff(p(x), x);
ddp := diff(p(x), x, x);
DELTA2 := piecewise(k <> j, -2*(-1)^(j-k)/(j-k)^2, k = j, -(1/3)*Pi^2)/h^2;
DELTA1 := piecewise(k <> j, (-1)^(j-k)/(j-k), k = j, 0)/h;
DELTA0 := piecewise(k <> j, 0, k = j, 1);
PHI := proc(x)
            ln(sinh(x));
       end proc;
dPHI := diff(PHI(x), x);
ddPHI := diff(PHI(x), x, x);
for i from -N-1 to N do
    x[i] := ln(exp(i*h)+(exp(2*i*h)+1)^(1/2));
 end do:

POL := seq(simplify(eval(sum(c[k]*((eval(2*dPHI*DELTA1), x = x[j])+eval(x[j]*ddPHI*DELTA1, x = x[j])+x[j]*(eval(dPHI^2, x = x[j]))*DELTA2), k = -N .. N)+eval(ddp, x = x[j])+2*(sum(c[k]*(eval(x[j]*dPHI*DELTA1, x = x[j])+DELTA0), k = -N .. N)+eval(dp, x = x[j]))/x[j]+(c[j]*x[j]+p(x[j]))^m, x = x[j])), j = -N-1 .. N):
 

dp := 2*x*c[-18]

ddp := 2*c[-18]

DELTA2 := 5.77847633136095*piecewise(k <> j, -2*(-1)^(j-k)/(j-k)^2, k = j, -(1/3)*Pi^2)

DELTA1 := 2.40384615384615*piecewise(k <> j, (-1)^(j-k)/(j-k), k = j, 0)

piecewise(k <> j, 0, k = j, 1)

proc (x) ln(sinh(x)) end proc

cosh(x)/sinh(x)

1-cosh(x)^2/sinh(x)^2

K := CodeTools:-Usage( fsolve({seq(POL[v] = 0, v = 1 .. 2*N+2)}) );

memory used=3.26GiB, alloc change=7.00MiB, cpu time=22.36s, real time=20.56s, gc time=4.10s

{c[-18] = -0.158067746087552e-1, c[-17] = -0.387542884266278e-4, c[-16] = -0.130630169518967e-3, c[-15] = -0.239834317398380e-3, c[-14] = -0.415811440194114e-3, c[-13] = -0.646693636830992e-3, c[-12] = -0.101389527576115e-2, c[-11] = -0.153641273761729e-2, c[-10] = -0.235331965565866e-2, c[-9] = -0.355711573437526e-2, c[-8] = -0.541208015126014e-2, c[-7] = -0.818412067244522e-2, c[-6] = -0.124116070087548e-1, c[-5] = -0.187405318724168e-1, c[-4] = -0.282444166168333e-1, c[-3] = -0.421629853858902e-1, c[-2] = -0.619121771983790e-1, c[-1] = -0.878051524903856e-1, c[0] = -.117654444653016, c[1] = -.145782460278171, c[2] = -.166307162557309, c[3] = -.176491199766487, c[4] = -.177378358017197, c[5] = -.171361063895020, c[6] = -.160969710441339, c[7] = -.147975310537005, c[8] = -.133706933631217, c[9] = -.118868260983916, c[10] = -.104020256360719, c[11] = -0.893402689793291e-1, c[12] = -0.750814767898076e-1, c[13] = -0.611879309401569e-1, c[14] = -0.478386383346360e-1, c[15] = -0.348219351557400e-1, c[16] = -0.224317459414568e-1, c[17] = -0.100575149091325e-1}

Digits,oldDigits := 50, Digits: # increase while computing residuals
eval(map(abs@(rhs-lhs),{seq(POL[v] = 0, v = 1 .. 2*N+2)}),K);
max(%);
Digits:=oldDigits: # restore

{0.1848364467526429241183563989099735e-15, 0.202587081085419172747222400227041933204473295e-15, 0.10097241581802646603811727298402707e-14, 0.14553137467375713934400384650335701e-14, 0.22560173770270976761107758930223125e-14, 0.27160112886445930168063402665723885e-14, 0.39319837003387081002098376695012516e-14, 0.51325880161426781211477890118403909e-14, 0.53776848377504449929291594039049726e-14, 0.585095713936770588566043075378528608079752e-14, 0.688270019244350727096392771883525748e-14, 0.78563785170572199376709531083519070e-14, 0.904099525662313748582996811317018268e-14, 0.92859555415356170929146484133354861e-14, 0.94853914917733280044366040147412311e-14, 0.97752062759665148221253544524779454e-14, 0.107763339630216080424904434349711800e-13, 0.11405979550333858018322263416696762e-13, 0.132198657572373080463771901878901716e-13, 0.145643280852386524219168072988804449e-13, 0.149782918036741699690104183753487886e-13, 0.152120601633383535853708830652785124275e-13, 0.163630089886033587689846874750524749e-13, 0.188390085778795556521433153435070662e-13, 0.225289266140715475831341813098011734e-13, 0.22643740696535725265433769519996728e-13, 0.24504005846030340232899949273655086e-13, 0.249295353634064239781347372712429555e-13, 0.2605477850150012807659458375862248648e-13, 0.268800255486011004914803275475734128e-13, 0.329272851091861856917522923708165545e-13, 0.3337556394391884128014195478415057e-13, 0.4857403341990820827906483739374796e-13, 0.578726405983914849751293908345002674e-13, 0.606234244112589733855448574102251748e-13, 0.7569341881763600803047512052987563131e-13}

0.7569341881763600803047512052987563131e-13

with(DirectSearch):

KDS:=CodeTools:-Usage( SolveEquations([seq(POL[v] = 0, v = 1 .. 2*N+2)], evaluationlimit = 10000000) );
 

memory used=13.61GiB, alloc change=0 bytes, cpu time=9.78m, real time=2.53m, gc time=52.76s

KDS := [2.90084071296831*10^(-8), Vector(4, {(1) = ` 1 .. 36 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}), [c[-18] = -0.158088096698655e-1, c[-17] = -0.387372604634768e-4, c[-16] = -0.130628354565356e-3, c[-15] = -0.239814730814512e-3, c[-14] = -0.415808213284412e-3, c[-13] = -0.646665453648750e-3, c[-12] = -0.101388696973837e-2, c[-11] = -0.153635877299765e-2, c[-10] = -0.235322045150403e-2, c[-9] = -0.355695760087773e-2, c[-8] = -0.541194604080468e-2, c[-7] = -0.818397616579811e-2, c[-6] = -0.124115320868477e-1, c[-5] = -0.187404175924928e-1, c[-4] = -0.282444401564308e-1, c[-3] = -0.421630930157411e-1, c[-2] = -0.619121766646340e-1, c[-1] = -0.878045849452816e-1, c[0] = -.117652147369853, c[1] = -.145779235574756, c[2] = -.166303042914011, c[3] = -.176488924141685, c[4] = -.177377976537771, c[5] = -.171360759204415, c[6] = -.160967434973583, c[7] = -.147971797900332, c[8] = -.133702754001262, c[9] = -.118863260180314, c[10] = -.104015284296128, c[11] = -0.893367200861269e-1, c[12] = -0.750793728225821e-1, c[13] = -0.611862816909835e-1, c[14] = -0.478375133919584e-1, c[15] = -0.348219354414091e-1, c[16] = -0.224325517898647e-1, c[17] = -0.100584010716369e-1], 381334]

Digits,oldDigits := 50, Digits: # increase while computing residuals
eval(map(abs@(rhs-lhs),{seq(POL[v] = 0, v = 1 .. 2*N+2)}),KDS[3]);
max(%);
Digits:=oldDigits: # restore

{0.1316051275014097210203166431394540633767326e-6, 0.16567493632487175641503270951110534929791808e-5, 0.171097061896518610704916484686124009437803975e-5, 0.198878038317539469674371030654542948189175e-5, 0.32964708225449765391612205407049708511591347820180e-5, 0.36325879158394784567603910731688966263904716e-5, 0.3851725808084988897920578176631049809733860e-5, 0.49646464490371853392283541801202920785065812e-5, 0.60890892790904263948088460983886903976605094e-5, 0.70363623448832614626242257806259121617999511e-5, 0.76145203964341754818177628290243078038844310242149e-5, 0.775220502813178515830217832601841447079846880e-5, 0.8378161538157416400386574710947981123823915e-5, 0.89421177524876502276344954543452592001914342e-5, 0.98170124856628187641607136453197188863566653e-5, 0.1169758246855752838450333464927366109118322984e-4, 0.139246505864902482964045769873769916140840802e-4, 0.140778123350927219393406017802690733781604448e-4, 0.166109274878519454023120195435751265679363930e-4, 0.169876715052390473458298370703925749175976160e-4, 0.178668817360814463379014807453011108554767105e-4, 0.181820800297383253512272665858361286004395892e-4, 0.188899101375935002010945680746218503476038793e-4, 0.188975806447715851676183514546955494103967398e-4, 0.2733275740365502174280626102023360593676048e-4, 0.276044871350656030202679618781482069944118799e-4, 0.280773598773883973037295194373159810278793796e-4, 0.28993793439522313932850202406749385764519814155465e-4, 0.327506996541199330313883762041458294139131200e-4, 0.330175317111639870334269695379056703090550518e-4, 0.3464188663803647672938358802772720534951015360e-4, 0.492112160463829717048186171468347841390459201e-4, 0.526231401385384954271656833023762174041184407e-4, 0.573621915174304647320863872247580496407557314e-4, 0.648075220779544423555542611235192153856057008e-4, 0.82731962706244942095321358821961365091919244e-4}

0.82731962706244942095321358821961365091919244e-4

 


 

Download bigSys_dig15.mw

 

 

Try adding the option,

tickmarks=[piticks,default]

to your call to implicitplot.

You might also wish to force the full horizontal view from 0 to Pi/2, say by also adding the option.

view=[0..Pi/2,default]

 

Execute the statement,

randomize():

at the start of your worksheet.

Q:=Matrix(Q, datatype=float[8]);

or,

Q:=Matrix(5,5,(i,j)->simplify(Q[i,j],zero));

or several other variants.

Similarly,

H:=Vector(H,datatype=float[8]);

or,

H:=Vector(5,(j)->simplify(H[j],zero));

or variants. As you've seen in other message threads, small nonzero imaginary components can be handled with `fnormal`, before calling `simplify`. Or you can be more forcing and use the `Re` command. 

nb.Without running your code in Maple, your Matrix K looks purely real, symmetric, and your Matrix M looks purely real, symmetric, positive-definite. If so then the eigen-solution should have all-zero imaginary components.

Are you using an older version of Maple? If I recall then in recent Maple version the zero imaginary components of complex[8] rtables are suppressed on printing.

The fnormal command let's you replace small floating-point values with 0.0 or 0.0+0.0*I, etc. This command allows you to control the fineness, optionally. See its Help page.

As Kitonum has shown, the command simplify(expr, zero) removes the 0.0 occurences (real or imaginary) from the expression `expr`.

 

Are you looking for plots:-implicitplot3d, where its first argument might be f(x,y,z)=K for K numeric?

Let's try an example with more corrner cases in it.

Below I use uneval quotes (single right-quotes) because neither typefunc nor last_name_eval is a protected name.

restart;

eq := diff(x(t), t) = t + x(t) - 1 + sin(x(t))
      + cos(t) + G(t, x(t), v(t)):

indets(eq,{name,'typefunc(name,Non(last_name_eval))'});

                               {t, v(t), x(t)}

indets(eq,'typefunc(name,Non(last_name_eval))');

                                {v(t), x(t)}

indets(eq, {name, function(name)});

                           {t, cos(t), v(t), x(t)}

indets(eq,function(name));

                            {cos(t), v(t), x(t)}

How about handling F(t), where F is a user-defined procedure which returns unevaluated if its argument is not numeric? (Think of using the known option of dsolve/numeric, and other similar situations.) For unassigned argument t, calling F(t) returns the unevaluated function call F(t) itself.

Since F is a known procedure, which is already defined to compute an output number from an input number, then neither F nor F(t) is a "dependent variable" of the DE.

# In general F may be a black-box procedure, whose workings
# are not explicitly known.
F:=proc(u)
  if not type(u,numeric) then
    return 'procname'(u);
  else
    tan(u)+u^2;
  end if;
end proc:

F(t); # returns unevaluated for unassigned t

                                    F(t)

F(0.1); # returns a number. we won't "solve for F".

                                0.1103346721

eq2 := diff(x(t), t) = t + x(t) - 1 +  sin(x(t))
       + cos(t) + G(t, x(t), v(t)) + F(t):

indets(eq2,{name,'typefunc(name,Non(last_name_eval))'});

                               {t, v(t), x(t)}

indets(eq2,'typefunc(name,Non(last_name_eval))');

                                {v(t), x(t)}

indets(eq2, {name, function(name)});

                        {t, F(t), cos(t), v(t), x(t)}

indets(eq2,function(name));

                         {F(t), cos(t), v(t), x(t)}

procedures (of which operators are a subset) are of type last_name_eval.

Try it with eval(f) instead of just f, in the line where you augment nonIdMaps.

Otherwise, as seen, you're just building up a list of repeats of the name `f`  (and by the time you get around to mapping f->f(y) over the list you've finished the loop, and all the f's evaluate to its last iterated value).

ps. `select` is a better choice for this task than is looping and augmenting a list, IMO. Say,

select(f->f(__y)<>__y,
       [x->x, x->2*x, x->3*x]);

I notice that all the code posted so far does not protect against y being assigned. Perhaps your situation is in a context where y is local and you know it is unassigned. I use __y just in case. You could also make y local to the predicate operator/procedure.

You could try DocumentTools:-Tabulate. The caveat is that command only allows for one such plot (or collection/Table of plots) printed immediately after each paragraph or execution group.

P1:=plot(sin(x),x=-Pi..Pi,size=[300,300]):
DocumentTools:-Tabulate([P1],alignment=left,exterior=none,
                        width=300,widthmode=pixels):
P2:=plot(sin(x),x=-Pi..Pi,size=[500,500]):
DocumentTools:-Tabulate([P2],alignment=left,exterior=none,
                        width=500,widthmode=pixels):

You can also do the above all manually, by inserting GUI Tables, toggling off their borders, making them left aligned, and using DocumentTools:-SetProperty to put plots into them programatically.

You can also force the output of a regular printing of a plot to be left-aligned using the main menubar. But as far as I know that quality is not persistent and is undone by re-execution. The same goes for plots:-display(Array([....]) where you weight the right portion with an empty plot and then toggle all the qualities of the inserted Table (which also don't persist upon re-execution).

By default the caret is too narrow, and is rendered to the left (perhaps not only because of the italic slant).

Two other ways to accomodate this are:

1) Augment the base symbol with an invisible space, so that the caret gets stretched and shifted right. This pushes the subscript slightly to the right, but the overall effect is still better than the original caret misplacement, IMO.

2) Use double-underscore rather than square-bracket indexing for the subscripting. This also shifts and strectches the caret. But (see below) the effect may be too much when the subscript is itself wide.

nb. In My Maple 2017.1 for Linux the carets looks a little narrower and sharper than they do as rendered inline below (by this site). But the placement and width effects are essential similiar.

# Doesn't look great, caret is too narrow and too left.
#
plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mi("sigma"),mo("&circ;"))`[y]]);

# More effort, but seems more robust (see below).
#
plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mrow(mi("sigma"),mi("&InvisibleTimes;")),mo("&circ;"))`[y]]);

# Looks ok here, but not in wider case (see below).
#
plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mi("sigma__y"),mo("&circ;"))`]);

#
# Now let's try those with a wider base as well as a wider
# subscript.
#

plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mi("&sigma;T"),mo("&circ;"))`[yyy]]);

plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mrow(mi("&sigma;T"),mi("&InvisibleTimes;")),mo("&circ;"))`[yyy]]);

plot([], size=[400,200], tickmarks=[[],[]],
     labels = [x, `#mover(mi("&sigma;T__yyy"),mo("&circ;"))`]);

 

Download overcaret.mw

Add the option phi = 0 .. infinity to your fsolve call, to avoid negative values for this example.

You can also use Student:-Calculus1:-Roots to get multiple roots within a finite range. Results may vary according to your Maple version, though. The first attachment below was run in Maple 2017.0, and you are using what, Maple 13.00 still?

ANALYTIC_1.mw

Here's an attachment that was run in Maple 13.00.

ANALYTIC_1_M13.00.mw

nb. Your original worksheet returned calls to Analytic, unevaluated, because you didn't first load the RootFinding package. You also failed to supply it with a (thin, complex) range for phi. [edited] See below for more, including another attachment by me.

Q2: For this example, you can obtain the desired simplification with just,

  simplify(ap2, size)

or,

  numer(ap2)/denom(ap2);

or,

  frontend(simplify, [ap2]);

Q3: note the t1(T) in the lhs of ap3 (as opposed to just your original t1 ).

ap3 := Diff(t1(T), T) = (m-T)/(m-t1(T));

ap4 := Diff(lhs(ap3), T) = subs(value(ap3), diff(rhs(ap3), T));

Q5: I don't understand why you thing that the originally supplied ap6 and ap7 are always equal.

Try the command `timelimit` which can wrap around another procedure call.

The `timelimit` command is intended for exactly the kind of situation you describe.

You can use the try...catch mechanism to control this flow.

When `timelimit` times out it emits an error, and your code can catch that, and proceed.

The command LinearAlgebra:-GenerateMatrix should be able to do that.

It's awkward to do, and not general, but... you can alter the Typesetting so that the exponent is smaller.

Making this more precise, say to affect only powers within sqrts, would be more involved.

Using labelfont=[...,bold,N] might make the two sizes appear more similar in weight.

restart;

N := 24:

plot(cos(1/theta)/sqrt(1-theta^2), theta=0.2 .. 0.99,
     size=[700,400],
     labels = ['theta', typeset(cos," (","1/",theta,") / ",
               sqrt(sort(1-theta^2,theta,ascending)))],
     labelfont=[Times,N],labeldirections=[horizontal,vertical]);

K:=eval('Typesetting:-Typeset'(sqrt(sort(1-theta^2,theta,ascending)))):
K:=subsindets(
  subsindets(K,specfunc({Typesetting:-mn,
                         Typesetting:-mo,Typesetting:-mi}),
  uu->op(0,uu)(op(uu),size=sprintf("%d",N))),
  specfunc(Typesetting:-msup),
             u->op(0,u)(op([1,0],u)(op([1,..],u),size=sprintf("%d",N)),
                        op([2,0],u)(op([2,..],u),size=sprintf("%d",ceil(0.5*N))))
  ):

plot(cos(1/theta)/sqrt(1-theta^2), theta=0.2 .. 0.99,
     size=[700,400],
     labels = ['theta', typeset(cos," (","1/",theta,") / ",K)],
     labelfont=[Times,N],
     labeldirections=[horizontal,vertical]);

 


cust_label.mw

First 194 195 196 197 198 199 200 Last Page 196 of 336