Carl Love

Carl Love

28110 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

You say that you increased the precision in the Options menu. What did you increase it to? If I increase it to 40, then I get accurate results (for the plot in question). Rather than using the options menu, it is better to use the Digits environment variable.

Digits:= 40:

You can test the accuracy of the computation by using the shake command. For example, let's apply it to the last point in your plot. The second argument to shake is the number of Digits of precision.

expr:= eval(m(450, (1/250)*x, 250), x= 100): #Make sure to use exact rationals in this line.
shake(expr, 10);

     INTERVAL(-3.44691114182*10^28 .. 3.45185884337*10^28)

shake(expr, 40);

    INTERVAL(1.96553109914429154721697440528449547786023 ..
    2.03451842008920754181567678008881546299414)

shake(expr, 50);

     INTERVAL(2.000000043280774985094379521387152201289707358740187 ..  
     2.000000043287679818113471795018001692671257692766963)

The above indicates that at 10 digits precision (the default), the computation is total garbage; at 40 digits precision, there are two accurate digits; and at 50 digits precision, there are 11 accurate digits.

In Maple's 2D input, it's illegal to put a space between a function's name and its arguments. You have, for example, printf ("*"). That needs to be changed to printf("*")

After you correct that, you still have a small logic error that prevents the first line from being printed. I'll let you figure that out. Let me know if you can't.

plots:-inequal(
     {cos(x) <= y and y <= sin(x) or sin(x) <= y and y <= cos(x)},
     x= 0..4*Pi, y= -1..1, xtickmarks= "piticks", scaling= constrained, color= red
);

nprintf("%5.3f", a);

It's possible, although a bit tricky, to program this so that the 5 and the 3 are not hardcoded. I'll work on it if you need that.

Dr Venkat Subramanian: Here's my Lobatto procedure. The amount by which I increase Digits (essentially, by n^(2/3)) may be excessive, but high-degree polynomial evaluations are notoriously prone to round-off error. Nevertheless, this gives accurate results with very few integrand evaluations.

 

restart:


Lobatto:= proc(f::algebraic, R::name=range(realcons), n::{posint, Not(identical(1))})
local
     x:= lhs(R), a, b,
     P:= (D@@(n-1))(x-> (x^2-1)^(n-1)/2^(n-1)/(n-1)!), DP:= D(P)(x),
     F, oldDigits, r
;
     oldDigits:= Digits;
     Digits:= Digits+1+ilog2(Digits)+iroot(n,3)^2;
     (a,b):= op(evalf(op(2, R)));      
     F:= unapply(eval((b-a)*f, x= (b-a)/2*x + (a+b)/2)/P(x)^2/n/(n-1), x);
     r:= (b-a)*(eval(f, x= a)+eval(f, x= b))/n/(n-1) + add(F(x), x= [fsolve(DP)]);
     evalf[oldDigits](r)
end proc:


Digits:= 15:


Lobatto(sin(sin(x)), x= 0..Pi, 14);

1.78648748195005

That's with 14 integrand evaluations. Let's count how many evaluations `evalf/Int` uses in procedure mode.

F:= proc(x) option remember; `if`(x::realcons, sin(sin(x)), 'procname'(args)) end proc:

evalf(Int(F, 0..Pi));

1.78648748195005

nops(op(op(4, eval(F))));

109

Conclusion: The Lobatto procedure might be worthwhile if the integrand is smooth and its evaluations are expensive.

 

Download Lobatto.mw

Your model is a linear function of the unknown parameters ab, and c, so you should use LinearFit. The fact that it's a nonlinear function of the independent variable t is irrelevant.

Statistics:-LinearFit([1, t, t^2], X, Y, t);

Nonetheless, NonlinearFit will still work. Unlike the other respondents, I was able to load and execute your entire worksheet (in Maple 16), including the NonlinearFit command, without error and without making any changes. I do however agree with the other respondents that you should never embed a large dataset like this in a worksheet. It should be read in from a file.

The result was

One way would be to use three colors, like this:

plot([seq([cos(t), sin(t), t= 2*Pi/3*(k-1)..2*Pi/3*k], k= 1..3)], color= [green, yellow, red]);

x=1;
map(`*`, %, 2);
map(`^`, %, 2);
or map2(`^`, 2, %); (It's not clear which you meant.)
map(`^``, %, 1/3);

All the above work the same way for inequalities (even when they are not valid for inequalities---be careful!).

The most important commands for manipulating equations and inequalities are rhs and lhs, which extract the right side and the left side, respectively.

The vast, vast majority (I'd guess more than 99%) of Maple procedures are not mentioned in the help. That doesn't mean that they are in the kernel. It usually means that they were not intended for use by end users. The procedure in question, TRDpolynomial_ring, is a local of module RegularChains. That means that it was only intended to be called from other code within the module. To see its code:

kernelopts(opaquemodules= false):  #Allow access to module locals.
showstat(RegularChains:-TRDpolynomial_ring);

restart:
alias(X= RootOf(z^8+z^6+z^5+z^3+1)):
int_to_poly:= proc(n::nonnegint, x::{name, specfunc(anything, RootOf)})
local k, B2:= convert(n, base, 2);
     add(B2[k]*x^(k-1), k= 1..nops(B2))
end proc:

A:= <140, 155, 162, 64;
     218, 12,  245, 50;
     36,  251, 34, 253;
     171, 251, 184, 37>:

B:= int_to_poly~(A, X);
 

For most uses of the new Matrix B, you'll need to append mod 2 to the command. For example, to compute the determinant,

Det(B) mod 2;

You can let plot choose the number of points and their positions. This is called adaptive plotting. For example, using the "more complex" oscillator that you gave in a Reply---the one that required numpoints= 100*tmax---do

U:= eval(u[1](t), res):
V:= eval(v[1](t), res):
plot([U(t), V(t), t= 0..tmax]);

You'll get essentially the same plot as you got with odeplot using numpoints= 100*tmax. The number of points chosen by plot is 6293. The problem with this command is that it takes a little more than an hour to produce the plot. Unlike odeplot, the t-values used by plot aren't evenly spaced.

plot([seq([i, sin(Pi/12*(i-1))], i= 1..12)], style= point, symbol= diamond, symbolsize= 9);

You need a percent sign with the T, like this:

a.a^%T;

A plus sign with no percent sign also works:

a.a^+;

You asked:

Does [Length of output exceeds limit of 1000000] mean that the function couldn't be solved and I have to modify the mathematical model?

No, it doesn't mean that. It means that the solution has been found, but is just too long to display.

If the results are just too long to be displayed, how can I see the results?

You should seriously reconsider your desire to see the results. They could be hundreds of pages long. I have seen one over 4000 pages long. Once you start the display, there is no way to stop it other than by killing the Maple process. The stop button does nothing in this case.

And, practically, what would you do with such a solution?

You should use definite integrals, not indefinite integrals that you then evaluate between limits of integration.

Examples:

int(int(r^2*cos(theta), theta= Pi/4..Pi/3), r= 15..18);

or

int(r^2*cos(theta), [theta= Pi/4..Pi/3, r= 15..18]);

First 239 240 241 242 243 244 245 Last Page 241 of 395