Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 362 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

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.

@Axel Vogt 

I don't understand the purpose of your reformulation. Yes, it is interesting that Maple can perform that symbolic integration of the MeijerG expression in an instant. Perhaps that was your point. But the Int/diff expression in Question doesn't use the integral of the MeijerG expression; rather, it uses the integral of the square of the derivative of the MeijerG expression. This is much more complicated.

@Axel Vogt 

Here is my simplest "clean" interpretation of the Int/diff expression. By "clean", I mean that each of the three nested definite integrals has its own dummy variable of integration which is not used outside that integral:

H3p_Clean:= proc(ms::algebraic, b::numeric, eps::positive:= 0.5*10^(1-Digits))
local x:= indets(ms, And(name, Not(constant))), D, X, Y;
     if nops(x)>1 then
          error "1st argument must be an expression of at most 1 variable"
     end if;
     x:= x[];
     D:= diff(ms,x);
     evalf(
          subsindets(
               Int(
                    diff(
                         subs(x= Y, D)*Int((-b*X+1)*Int(D^2, x= 0..X), X= 1..Y),
                         Y
                    )*subs(x= Y, ms),
                    Y= 0..1
               ),
               specfunc(anything, Int),
               J-> Int(op(J), epsilon= eps)
          )
     )
end proc:

(I hope that my indentation, line break, and parenthesis style serves to elucidate rather than obfuscate the expression for you.)

In all cases that I've tried, this produces the same numeric results as the original "sloppy" version (posted in the Answer above). (By "sloppy", I mean that the same variable is used for all integrations and differentiations, and it's also used in the limits of the inner integrals.) But the weird thing is that in all cases the sloppy version ran significantly faster, usually by a factor of 2. So maybe the sloppy version is better.

I tried all of the following for the first parameter, msxx^2x^3sin(x)cos(x)exp(x), exp(-x^2)sin(x)+cos(2*x). For each, I used both and -1 for the second parameter, b. I also used the original MeijerG expression with b = .7. I varied the order in which the two integrals were done.

I have Digits set to 15, but I never used an epsilon smaller than .5e-10.

You could do an animation where each frame was one of a sequence of parallel planes intersecting the solid. In each plane, color would represent the value of the function. If you've ever seen an animation of a CAT scan or MRI, you'll know what I mean.

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