Carl Love

Carl Love

28070 Reputation

25 Badges

13 years, 36 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@MDD You simply cut and paste and execute that code that Roman gave, then 

with(PolynomialIdeals):
EliminationIdeal(<x>,{a,b,c});

@Markiyan Hirnyk 

By what criterion are you saying that the result produced by DirectSearch:-DataFit is better than that produced by my RecenteredPowerFit (which calls Statistics:-PowerFit)? You can't directly compare the residuals as reported by the various fitmethods because they are each computed by their own formula. If you compare the sum of the squares of the residuals, my fit is several times better than yours. It is also clear from the plots that my fit is better than yours. (My fit is the red curve in my plot; ignore the blue curve.)

Any result with > 1 is completely ridiculous, not just "not so good". Try plotting them. Of the nine fitmethods, seven produce such results for me (on the first run).

Furthermore, the results of DirectSearch:-DataFit are not entirely reproducible. Rerunning the exact same call (of course, without restart) often gives a completely different result, often changing a decent result into a completely ridiculuous one.

The vote up may be from the OP and may be due to the fact that I worked back and forth with her over three days and that I produced code that she was able to execute and with which she was able to get what she considered to be "perfect" results.

 

@Jerome BENOIT I also don't care what the shown material is called. I just didn't want other readers to get the impression that it was ever possible to see the complete code of a module---not by eval or by any other means. Your usage of "body" was likely to give that false impression.

Here's a suggestion: stabilize the plot to the unit circle by dividing by v; get rid of the view option. This makes animations that vary n or m much easier. So here's my version, a modification of your cycler from Reply "without t", not the one with the circles. I also added _rest to the plot command to make it easier to pass in plot options.

cycler := proc(k, p, m, n, T) local expr, t, u, up, v;
     u := exp(k*I*t);
     up := exp(-k*I*t*p);
     expr := exp(I*t) * (1 - u/m + I*up/n);
     v := 1 + abs(1/n) + abs(1/m);
     plots:-complexplot( expr/v, t = 0 .. T, axes = none, _rest);
end proc:

And here's an animation varying n. Note that I vary it exponentially.

plots:-animate(
     cycler, [5, 3, 2, 10^t, 2*Pi, color= purple], t= -2..1.5,
     frames= 100, paraminfo= false
);

And here's one varying m, also exponentially:

plots:-animate(
     cycler, [5, 3, 10^t, 3, 2*Pi, color= orange], t= -2..2.5,
     frames= 100, paraminfo= false
);

@Jerome BENOIT 

No version of Maple has ever shown a module's body, which is the executable code after the declarations. Indeed, that code is discarded after it's executed; it's not stored in the module. What eval shows (or used to show) is the module's header or shell, which includes the description, the options, and the names (but not the values, if any) of the locals and exports.

@Alejandro Jakubi What if you set interface(verboseproc= 2)?

@MDD 

module is a Maple data structure similar to a procedure that allows for a special type of local variable called an export which can be accessed outside the module. It is the primary structure in Maple for object-oriented programming, similar to a class in C++ or Java. Often, the exports of a module are procedures which allow controlled access to the locals. See ?module. Sometimes a module is used primarily or even exclusively as a container object like a table. See ?record.

This has nothing to do with the algebraic structure known as a module.

@Haley_VSRC 

Markiyan is using the DirectSearch package, which must be downloaded from the Maple Applications Center. It's a very good package with many different algorithms for nonlinear fitting. But, as you can see, several of those algorithms return "crazy" results with > 1.

@Haley_VSRC 

Here's how to use my procedure. In the final graph, note that there's a trend in the rightmost quarter or so of the data that is not explained by either curve.


restart:

TW:= <<2.00E-05|0.00723>,<4.86E-05|0.033171>,<7.31E-05|0.047349>,<9.77E-05|0.057941>,<0.000122|0.066835>,<0.000147|0.074711>,<0.000171|0.081864>,<0.000196|0.088463>,<0.000221|0.094621>,<0.000245|0.10042>,<0.00027|0.10591>,<0.000295|0.11114>,<0.000319|0.11614>,<0.000344|0.12091>,<0.000368|0.12544>,<0.000393|0.12977>,<0.000418|0.13391>,<0.000442|0.13788>,<0.000467|0.1417>,<0.000492|0.14539>,<0.000516|0.14896>,<0.000541|0.15244>,<0.000565|0.15582>,<0.00059|0.15911>,<0.000615|0.16231>,<0.000639|0.16544>,<0.000664|0.1685>,<0.000689|0.17148>,<0.000713|0.1744>,<0.000738|0.17725>,<0.000762|0.18005>,<0.000787|0.18279>,<0.000812|0.18548>,<0.000836|0.18811>,<0.000861|0.1907>,<0.000885|0.19324>,<0.00091|0.19574>,<0.000935|0.19819>,<0.000959|0.2006>,<0.000984|0.20297>,<0.001008|0.20531>,<0.001033|0.20761>,<0.001058|0.20987>,<0.001082|0.2121>,<0.001107|0.2143>,<0.001131|0.21647>,<0.001156|0.2186>,<0.00118|0.22071>,<0.001205|0.22279>,<0.00123|0.22484>,<0.001254|0.22687>,<0.001279|0.22886>,<0.001303|0.23083>,<0.001328|0.23278>,<0.001352|0.2347>,<0.001377|0.2366>,<0.001401|0.23847>,<0.001425|0.24032>,<0.00145|0.24215>,<0.001474|0.24396>,<0.001498|0.24575>,<0.001523|0.24751>,<0.001547|0.24926>,<0.001571|0.25099>,<0.001596|0.25277>,<0.00162|0.25433>,<0.001644|0.25571>,<0.001668|0.25708>,<0.001693|0.25843>,<0.001717|0.25978>,<0.001741|0.26112>,<0.001765|0.26245>,<0.001789|0.26377>,<0.001813|0.26508>,<0.001838|0.26638>,<0.001862|0.26767>,<0.001886|0.26895>,<0.00191|0.27022>,<0.001934|0.27148>,<0.001958|0.27273>,<0.001982|0.27396>,<0.002006|0.27519>,<0.00203|0.27641>,<0.002054|0.27762>,<0.002078|0.27881>,<0.002102|0.28>,<0.002126|0.28118>,<0.00215|0.28235>,<0.002174|0.28352>,<0.002198|0.28467>,<0.002222|0.28582>,<0.002246|0.28695>,<0.00227|0.28808>,<0.002294|0.2892>,<0.002318|0.29032>,<0.002342|0.29142>,<0.002366|0.29252>,<0.00239|0.29361>,<0.002414|0.29469>,<0.002438|0.29577>,<0.002462|0.29684>,<0.002486|0.2979>,<0.002509|0.29895>,<0.002533|0.3>,<0.002557|0.30104>,<0.002581|0.30207>,<0.002605|0.3031>,<0.002629|0.30412>,<0.002653|0.30513>,<0.002677|0.30614>,<0.002701|0.30714>,<0.002724|0.30813>,<0.002748|0.30912>,<0.002772|0.3101>,<0.002796|0.31107>,<0.00282|0.31204>,<0.002844|0.313>,<0.002868|0.31396>,<0.002891|0.31491>,<0.002915|0.31586>>:

RecenteredPowerFit:= proc(TW::Matrix, b_range::range(realcons))

local

     T:= TW[..,1], W:= TW[..,2],

     b:= Optimization:-NLPSolve(

          b-> Statistics:-PowerFit(T-~b, W, output= residualsumofsquares),

          b_range, method= branchandbound, objectivetarget= 0

     )[2][1],

     Res:= Statistics:-PowerFit(T-~b, W)

;

     ['A' = Res[1], ':-b' = b, 'c'= Res[2]]

end proc:

 

First try with the full data set. Smallest t value is 2e-5.

Sol_Carl1:= RecenteredPowerFit(TW, -1..TW[1,1] - 1e-7);

[A = HFloat(5.028381968959027), b = HFloat(1.9211938376704556e-5), c = HFloat(0.46614083419427993)]

Sol_Haley:= [A = 2.48, b = 1.1e-4, c = .35]:

PowerFitCheck:= proc(TW::Matrix, b::realcons)
local
     T:= TW[..,1] -~ b, W:= TW[..,2], n:= op([1,1],TW),
     nRSS:= Statistics:-PowerFit(T, W, output= residualsumofsquares)/n
;
     nRSS, Statistics:-PowerFit(T, W)
end proc:
     

PowerFitCheck(TW, eval(b, Sol_Carl1));

0.844786726509007e-3, Vector(2, {(1) = 5.02838196895903, (2) = .466140834194280})

Res:= PowerFitCheck(TW[5.., ..], eval(b, Sol_Haley));

Res := 0.215501674137240e-2, Vector(2, {(1) = 2.49659481877520, (2) = .354232604162157})

Sol_Haley5:= [A= Res[2][1], Sol_Haley[2], c=Res[2][2]];

[A = HFloat(2.496594818775202), b = 0.11e-3, c = HFloat(0.3542326041621568)]

Sol_Carl5:= RecenteredPowerFit(TW[5.., ..], -1..TW[5,1] - 1e-7);

[A = HFloat(3.682236898705571), b = HFloat(6.28308061404132e-5), c = HFloat(0.4158452235037357)]

PowerFitCheck(TW[5.., ..], eval(b, Sol_Carl5));

0.139855979031156e-3, Vector(2, {(1) = 3.68217591144777, (2) = .415842527560050})

P1:= plot(TW, style= point, symbol= cross, legend= data):

model:= A*(t-b)^c:

P:= plot([eval(model, Sol_Carl5), eval(model, Sol_Haley5)], t= 0..0.0029, legend= [Carl, Haley]):

plots:-display(P1,P);

 


Download RecenteredPowerFit.mw

@tomleslie The procedure that throws the error, `dsolve/numeric/DAE/make_proc`, is a monster: 792 lines and 36 parameters (as of Maple 18). The problem starts when floats in the equation are converted to exact rationals in line 22 and ends with the error message in line 719. I haven't traced what happens in between---it's a nightmare. It looks like the procedure has become a patchwork quilt over the years as features have been added to dsolve(..., numeric)

The same bug is manifested if the coefficient 3 is changed to 3.0.

@Haley_VSRC Would you please check the value of c that you transcribed above? I think that you may be off by a decimal place. I concur on your value of (given your b and treating the first four points as outliers), but I get c = 0.35.

Anyway, > 1 is impossible, since the data are concave down.

If you use my method and make just the first point an outlier, you'll get

[A = 4.32368499833578, b = 0.335082182570315e-4, c = .442173738481387].

Using my method decreases the average squared residual by a factor of more than 6 over your method.

Very nice. I love ornate pentagrams. I've figured out what the parameters and do just by inspection. What do pm, and n do? And why is t a parameter rather than a local? Surely it makes no sense for t to be anything other than a name.

@MDD I use the functional operator form (the arrow form) when a procedure is simple enough that it qualifies for it. The procedure Variables above is essentially equivalent to 

Variables:= proc(e, p::{set,list,rtable}(name)) 
     indets(e, And(name, Not(constant))) minus convert(p,set)
end proc:

just like f:= x-> x^2 is essentially equivalent to f:= proc(x) x^2 end proc.

See help page ?operators,functional.

@lham You wrote:

for w^0 , is it correct ?
coeff(ex,w,0)/w^0+coeff(ex, w, -1)/w + coeff(ex, w, -2)/w^2;

Yes, it is correct. The expression in Question doesn't have any w^0 terms (also called constant terms), but, if it did, that would be the way to extract them.

-1 or -2 in above what it means ?

It's the exponent of (since 1/w^2 = w^(-2)).

and another question , i want to write original equation without any factoring 

Use the expand command.

 

Please put your Questions in the Questions section, not the Posts section. I moved this one for you.

First 483 484 485 486 487 488 489 Last Page 485 of 709