Carl Love

26832 Reputation

11 years, 287 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

Make each frame a seq of all previous fr...

The trick is to make each frame include a static (i.e. not insequence) display of all previous frames. This can be done with any animation. Here's the original code (which I indented) and my modifications. I added a splash of color. I also changed the sum to add for efficiency, which necessitated making g a procedure.

(***** begin original (comment out) *****
g:= sum(arctan(1/sqrt(n)), n = 1 .. i):
p:= seq(
plots[polygonplot](
[[0, 0],
[sqrt(k+1)*cos(subs(i = k, g)), sqrt(k+1)*sin(subs(i = k, g))],
[sqrt(k)*cos(subs(i = k-1, g)), sqrt(k)*sin(subs(i = k-1, g))]
], color = white
), k = 1 .. 70
):
plots[display]([p], insequence = true);
***** end original (comment out) *****)

N:= 70:
p:= [seq](
plots:-polygonplot(
[[0, 0], sqrt(k+1)*~[cos,sin](g(k)), sqrt(k)*~[cos,sin](g(k-1))
], color= COLOUR(HUE, k/N)
), k= 1..N
):
plots:-display([seq](plots:-display(p[1..K]), K= 1..N), insequence = true);

How to use evalhf...

The trick to using evalhf here is that the integration needs to be done at the time h is created, not at the time that h is called. And there's extra complexity that n = 2 needs to be handled as a separate case to avoid division by zero. So, you need to define h this way:

h := subs(J= int(f(x)*sin(n*x), x= -Pi..Pi)/Pi, n-> `if`(n=2, 1, J));

Your f and g can remain the same. If you define h as above, then the call evalhf(g(2)) will work. This is correct usage, not effective usage. To use evalhf effectively, you need to pack as much computation as possible into each invocation of evalhf.

Make the order 2...

I can't answer the "why?" of your question, but there is a very easy workaround: Make the order 2 or more.

No "freeze" for integers and rationals....

Unfortunately, freeze doesn't work for integers and rationals. But subs does work. So the following kludge does it:

ex:= 16^j/16;
exf:= subs(_x= 16, simplify(subs([1/16= 1/_x, 16= _x], ex)));

It sure looks ugly though.

Extra asterisk...

You have a line

H := dsolve*(ivp, {e(t), w(t)}, numeric);

The asterisk in that line should be removed.

Also, it is not an error, but there is no good reason to use `union` in prefix form. You might as well write

ivp:= sys union ic;

which is closer to the usual mathematical notation.

There are many ways to plot it. Here are two:

plots:-odeplot(H, [[t,e(t)],[t,w(t)]], t= 0..100, color= [blue,red]);

plots:-odeplot(H, [e(t),w(t)], t= 0..100);

Why an array of functions?...

Why do you want it to be an array of functions? Why not just have phi be a function of two variables, i and x?

phi:= (i,x)-> (i+2)/(i+1)*x^i + x^(i+1);

If you want to get fancy with the syntax and still be able to invoke phi with the square brackets, you can do this:

phi:= proc(x)
local i;
i:= op(eval(procname));
(i+2)/(i+1)*x^i+x^(i+1)
end proc:

phi[2](x);

Both of these approaches have the advantage of not polluting the memory with n copies of what is essentially the same function.

What do you mean by the color .1?...

This works:

DE := [diff(x(t), t) = sin(x(t))/sin(t)];
with(DEtools):
DC:= sin(x(t))/sin(t):
dfieldplot(DE, x(t), t = -2 .. 2, x = -1 .. 2, arrows = medium, title = 'Diff', color= DC);

Note that x is always x(t) in the differential equation. What were you trying to achieve with the .1 in the color?

solve, h?, not so simple...

The command is solve(..., B). It is not clear what you mean by h and by "equation". Is h meant to be a variable in the equation? Then your := should be =. Or do you mean that the expression assigned to h should be set to 0? Either way, the answer is not so simple: You have a degree-6 polynomial in B.

Split cases via |A intersect B| = 3, 4, ...

Notation: I'll use juxtaposition to indicate intersection, `+` for union, `-` for minus, with `+` having precedence. |A| is nops(A), and (n C k) is a binomial coefficient.

We split cases via |AB|. There are three cases, |AB| = 3, 4, or 5. Regardless of the case, we have a factor of (30 C 5) for the choice of A+B and a factor of 2^25 choices for C - A+B.

Case |AB| = 5: Thus A=B. There are (5 C 3) = 10 choices for ABC in AB.

Case |AB| = 4: Either |A|=4 and |B|=5 or vice versa. This dichotomy contributes a factor of 2, and we will consider only the |A|=4 case without loss of generality. There are (5 C 4) choices of A within A+B and (4 C 3) choices of ABC in A. Total factor for this case 2*(5 C 4)*(4 C 3) = 40.

Case |AB| = 3:

• Subcase |A|=3, |B|=5 or vice versa. Factor of 2 for the dichotomy. (5 C 3) choices of A in A+B, and only 1 choice of ABC in A. Total factor: 20
• Subcase |A|=4, |B|=4. (5 C 4) choices of A in A+B, (4 C 3) choices of AB in A, and(3 C 3) choices of ABC in AB, for a total of 20.

Grand total: (20+20+40+10)*(30 C 5)*2^25 = 430353709793280 =

Edit: Corrections to the above.

Case |AB| = 4: Either |A|=4 and |B|=5 or vice versa. This dichotomy contributes a factor of 2, and we will consider only the |A|=4 case without loss of generality. There are (5 C 4) choices of A within A+B and (4 C 3) choices of ABC in A. The choice of whether the lone element of B-A is or is not in C contributes a factor of 2. Total factor for this case 2*(5 C 4)*(4 C 3)*2 = 80.

Case |AB| = 3:

• Subcase |A|=3, |B|=5 or vice versa. Factor of 2 for the dichotomy. (5 C 3) choices of A in A+B, and only 1 choice of ABC in A. For each of the 2 elements of B-A, there is a choice of whether it is in C, for a factor of 2^2. Total factor: 2*(5 C 3)*2^2 = 80.
• Subcase |A|=4, |B|=4. (5 C 4) choices of A in A+B, (4 C 3) choices of AB in A, and(3 C 3) choices of ABC in AB. For each of the 2 elements of (A-B)+(B-A), there is the choice of whether it is in C. Total factor: (5 C 4)*(4 C 3)*2^2 = 80.

Grand total: (10+80+80+80)*(30 C 5)*2^25 = =

floating point contagion...

There must have been a floating-point number introduced into the computation, either explicitly or through use of evalf or a similar command. What were the steps that led to your 2.*I*q*r? This is a phenomenon called "floating-point contagion", (which makes it sound like a bug, but it it the intended behavior). Basically, the decimal points spread from number to number in your expression. To get rid of them,

subsindets(ex, float, convert, rational);

where ex is the expression that you want to correct.

algsubs, not subs...

You can't use subs to substitute for something that doesn't exist as a distinct entity. Notice the difference in the following two uses of subs.

 > restart;
 > ex:= (x+y)^2 + (a+b+c)^2;

 > subs(x+y= 2, ex);

 > subs(a+b= 2, ex);

To be precise about what I mean by "distinct entity" requires some knowledge of Maple's internal representations of objects.

But, we can ignore that, because there is a command that is much more sophisticated than subs, and which doesn't care about the internal representations. That command is algsubs. It only does one substitution at a time. So, to use it for your case, we need to put it into a loop.

Replace the line

Eqs1:= subs(BC, Eqs);

with

for s in BC do Eqs:= algsubs(s, Eqs) end do; Eqs1:= Eqs;

I did this in your code, and it ran through to the end, generating the point plot.

sort, numer, denom...

Assuming that ex is set to the expression that you want to manipulate, it can be done as

sort(-numer(-ex)/denom(ex), x);

The second argument to sort, the x, is to tell it which variable that you want to appear first.

select(has, G1, w)...

I can show you how to do it by TypeTools if that's what you really want. But isn't this preferable?

select(has, G1, w);

Let me know.

Also, you are using the wrong quotes again, you mean 'w', not `w`. The latter does nothing to something which is already a name. As for the '[]', I doubt that that's what you really want. There's `[]`, [], and [][]. They all do different things when applied as functions.

By integration with a correction factor...

I got 5/4*5^(m+1)/m, but that is asymptotically equivalent to your result (since the limit at infinity of their ratios is 1). I did it by integrating, getting the series for that, and then applying a correction factor.

 > restart;
 > interface(typesetting= standard): interface(prettyprint= 3):
 > Order:= 2:
 > assume(m > 0);
 > s:= j-> 5^j/j:

Verify that getting the asypmtotic series the normal way won't work.

 > S1:= sum(s(j), j= 1..m+1);

 > asympt(S1, m);

Error, (in asympt) unable to compute series

The summand s(j) is increasing for j > 1. Imagine approximating the integral of s(j) with rectangles under the curve. Then the sum S1 would represent the integral on 1..m+2. Likewise, the integral can approximate the sum.

 > int(s(j), j= 1..m+2);

 > asympt(%, m);

 > A1:= op(1,expand(%));

Estimate the error factor: that is, the factor by which the integral over the last rectangle overestimates that rectangle.

 > int(s(x)/s(m+1), x= m+1..m+2);

 > asympt(%,m);

Since this is asymptotic to O(1), we are justified in applying a constant correction factor A1. So, the final first term of the asymptotic series for S1 is

 > A:= A1/op(1,%);

Verify with a plot.

 > plot(S1/A, m= 30..1000);

Some CDF identities...

First, I'll create Z and its CDF again. (CDF stands for Cumulative Distribution Function, if that helps.) But, you only need to do this once per session. This Z is the same N(0,1) as in your first problem.

Z:= Statistics:-Distribution(Normal(0,1)):

And I'll also define Φ for convenience (less to type):

Phi:= evalf@Z:-CDF:

You can use the above for all your "Z" problems.

Now, to the problem at hand. First some identities that I hope that you already know. If you don't know these, find them in your textbook and memorize them.

1. For any random variable (r.v.) Z and any real z ≥ 0, P(|Z| > z) = P(Z > z) + P(Z < -z).
2. For any continuous r.v. Z and any z, P(Z ≥ z) = P(Z > z), (i.e. isolated points don't matter).
3. For any continuous r.v. Z and any z, P(Z > z) = 1 - P(Z < z).
4. For Z ~ N(0,1) and any z, P(Z > z) = P(Z < -z) (symmetry).