## 7776 Reputation

17 years, 361 days

## Animate pursuit...

That's a rather challenging problem.  It took me some time to figure it out, and therefore I wrote down the details so that I won't forget.

But...

When I display my worksheet here, it shows a lot of garbage.  I don't know why.  You may want to download the worksheet from the link below and view it in Maple.  This anumation is extracted from the worksheet.

## No solution?...

It appears that your system of equations has no solution. To see that, plot zz := rhs(eqE)-lhs(EqE) as a function of E1 and E2, with a variety of choices for nu and h1 within the specified ranges.  You will see that the zz is always negative, meaning that eqE cannot be satisfied.

```zz := unapply((rhs-lhs)(eqE), nu, h__1);

plot3d(zz(0.1, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 500), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.1, 1000), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.3, 1000), E__1=0..5000, E__2=0..5000);

plot3d(zz(0.5, 1000), E__1=0..5000, E__2=0..5000);

```

## See if this does what you want. ...

See if this does what you want.

 > restart;
 > w:=(1-delta[h]-delta[b])*(1-phi/(kc-3/2))^(-kc+3/2)*(-kc+3/2)/((kc-3/2)*(1-phi/(kc-3/2)))+delta[h]*(1-phi/(sigma[h]*(kh-3/2)))^(-kh+3/2)*(-kh+3/2)/((kh-3/2)*(1-phi/(sigma[h]*(kh-3/2))))-(1/18)*delta[b]*sqrt(3)*(3*sqrt((M-u+sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu-3*sqrt((M-u-sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu)/(mu*sqrt(mu*sigma[b]))-(1/18)*sqrt(3)*(-3*sqrt((M+sqrt(3)*sqrt(sigma[i]))^2-2*phi)+3*sqrt((M-sqrt(3)*sqrt(sigma[i]))^2-2*phi))/sqrt(sigma[i]);

 > e1 := eval([w,diff(w,phi)],phi=0);

 > e2 := simplify(e1) assuming positive;

 > e3 := simplify(e2) assuming positive,         M - u - sqrt(3)*sqrt(mu)*sqrt(sigma[b]) > 0,         M - u + sqrt(3)*sqrt(mu)*sqrt(sigma[b]) > 0,         M - sqrt(3)*sqrt(sigma[i]) > 0;

## Built-in...

Maple's OrthogonalExpansions library can calculate the Fourier series for you.

 > restart;
 > with(OrthogonalExpansions):
 > f := x -> piecewise(x<0, 0, x);

 > F := FourierSeries(f(x), x=-Pi..Pi, n);

 > plots:-animate(plot, [F, x=-Pi..3*Pi, color=blue], n=1..10, frames=10);

## A solution...

The answer to your Question #1 is easy.  Here it is. I haven't worked out the details of #2.  Perhaps someone else will.

 > restart;
 > with(plots):

Cylinder of radius  centered on the  axis; parameters  and

 > C1 := ;

Cylinder of radius centered on the axis; parameters  and

 > C2 := ;

Intersection curve

 > C1 =~ C2; solve(%, {s,p,q}); curve := eval(C2, %);  # parametrized by t

 > %display(         %tubeplot([seq](curve), t=-Pi..Pi, radius=0.03, color=yellow),         %plot3d(C1, s=0..Pi, p=-a..a, color=blue),         %plot3d(C2, t=-Pi..0, q=0..1.5, color=red, transparency=0.5),         %plot3d(sqrt(a^2-y^2)+1e-3, x = -sqrt(c^2-y^2) .. sqrt(c^2-y^2), y=-c..c, color=green), scaling=constrained, style=surface, axes=none, orientation=[20,65,0]): value(subs(a=1.1, c=1, %));

 >

## Elliptic section...

This worksheet shows how to calculate the parametric equations of the cone, the plane, and the ellipse that lies in their intersection.  The cases of parabola and hyperbola can be done in the same way.

Intersection of a plane and a cone

 > restart;
 > with(plots):

Parameters:

: half-angle of the cone's vertex

:  angle of the plane's normal relative to the  axis   (set  to get an ellipse)
:   coordinate of the plane's intersection with the  axis

 > params := beta=0.9*alpha, alpha=Pi/6, d=1/2;

Parametric equation of the cone ( and  are the parameters)

 > C := s*< tan(alpha)*cos(t), tan(alpha)*sin(t), 1>;

 > p1 := plot3d(subs(params, C), s=-1..1, t=-Pi..Pi,         scaling=constrained, color="Red", transparency=0.1);

The vectors  form an orthonormal triad, with  normal to the
plane and  within the plane:

 > N :=  < sin(beta), 0, cos(beta) >; A := < 0, 1, 0 >; B := LinearAlgebra:-CrossProduct(N,A);

Parametric equation of the plane. ( and  are the parameters)

 > P := A*u + B*v + d*<0,0,1>;

 > p2 := plot3d(subs(params, P), u=-1..1, v=-1..1, color="Blue", transparency=0.7);

To find the parametric equation of the ellipse , we equate the parametric equations

of the cone and the plane.  That gives us a system of three equations in the four unknowns

.  We solve the system for  in terms of : and then plug the result into the

equation of the cone to get the parametric equation of the ellipse, parametrized with :

 > C -~ P: seq(%): solve({%}, {u,v,s}): E := subs(%, C);

 > p3 := tubeplot([seq](subs(params, E)), t=-Pi..Pi, radius=0.02, color="Yellow");

 > display(p1,p2,p3,style=surface,scaling=constrained, axes=none);

## Try this simpler one...

If you don't know how so solve a problem, try a simpler one first.  See if you understand the following before you tackle your messy case.

 > restart;
 > with(LinearAlgebra):
 > Digits := 30;
 (1)

Make a random matrix; it is very likely to have a nonzero determinant

and therefore be nonsingular:

 > M__exact := RandomMatrix(4); Determinant(M__exact);
 (2)

Replace the (1,2)  element by , and determine  to make  into a singular matrix:

 > M__exact[1,2] := x: Determinant(M__exact): solve(%, x): M__exact[1,2] := %;
 (3)

Now  is singular:

 > Determinant(M__exact);
 (4)

Reduce the accuracy by truncating the entries from 30 to 20 decimal points:

 > M := evalf[20](M__exact);
 (5)

Now  is approximately singular,

 > Determinant(M);
 (6)

but from the point of view of Digits = 30,  is not singular, and therefore

the null space is empty:

 > NullSpace(M);
 (7)

and that's the issue that you are having in your worksheet.
To calculate the null space of the
approximately singular matrix ,

apply the singular value decomposition:

 > S, Vt := SingularValues(M, output=['S','Vt']);
 (8)

The returned vector  is the list of the singular values of .  We wee that

the last singular value is of the order of , which is almost

zero.  That says that the last row of the matrix  is the basis

of the corresponding approximate null space.  Let's check:

 > null__vec := Vt[-1]^+;
 (9)
 > M . null__vec;
 (10)

We see that the product is essentially the zero vector, which confirms

that  is in the null space.

## Corrections...

 > restart;

You are solving this initial value problem numerically:

 > f := r -> 1 - 1/r: eqq := {RR(2) = 2 + log(1), diff(RR(r), r) = 1/f(r)};

and then you write a (very inefficient) procedure for finding the inverse function
of that solution numerically.

Why not solve for the inverse right from the beginning?  It can be done
symbolically and it's very fast!  To explain, let me change your  notation

to a single-letter symbol . Then the differential equation is
that is .   Your initial condition is

 > de := diff(r(q),q) = 1 - 1/r(q);  ic := r(2)=2;

 > dsol := dsolve({de, ic});

There you have it.  That's equivalent to your numerical procedure

I will call it .

 > RT := unapply(rhs(dsol), q);

Plot your when you get the chance, and compare it with:

 > plot(RT(R), R=-4..12);

The rest of your worksheet works correctly as is.

## The proper demo...

Your attempt to animate the tan function misses an essential ingredient.

Consider the trigonometric circle (a unit circle centered at the origin in Cartesian coordinates) and a ray connecting the origin O = (0,0) to a point P = (cos(t), sin(t)).  Then, the horizontal projection of P is (cos(t),0) and the vertical projection is (0,sin(t)).  This enables us to associate the sine and cosine functions with the diagram's geometry.

Question: Where does the tan function fit in that diagram?

Answer:  Draw a tangent line L to the circle at (1,0).  Extend the ray OP to intersect the line L at some point H. You will see right away that H = (1, tan(t)).  So the tan function is associated with the point H in our diagram.  This is an absolutely essential information that should be conveyed to the students when teaching them the trigonometric circle.  [The function cot also has a similar geometric interpretation; I will leave it to you to figure it out.]

The animation below demonstrates the tan function properly in the context of the trigonometric circle.  Note the vertical tangent line to the unit circle.

AnimationCercleTrigoFonctiTrigo-rr.mw

## There is no neat solution...

 > restart; q1:=1-1/(w**2-3*sigma*k**2)-deltab*mu/((w-k*u0b)**2-3*mu*deltab*sigmab*k**2)+A/k**2=0;

You want to solve for w/k, so let .  Then , and the equation changes to

 > q2 := subs(w=Q*k, q1);

You want to solve eq2 for .  Let's try a very special case where most of the parameters are taken to be 1:

 > q3 := subs(deltab=1, sigmab=1, u0b=1, sigma=1, A=1, q2);

 > solve(q3, Q);

We see that even in this very special case, the solution  is a root of a fourth degree polynomial.

We conclude that there is no hope of obtaining a useful symbolic solution for the general case.

## A clean way...

Acer has offered an answer by modifying your worksheet as you have requested.  But your worksheet is much more complicated than it needs be.  Here is a clean way of producing your diagram.  You should be able to merge this with the improvements suggested by acer.

 > restart;
 > with(plots):
 > with(plottools):
 > angles := { seq(i*Pi/4, i=0..7), seq(i*Pi/6, i=0..11) };

Distance of angle lables from the center.  Adjust to taste.

 > r := 1.12;

 > display(         circle(color=blue, thickness=2),         seq(textplot([r,a,a],coords=polar), a in angles),         seq(line([0,0], [cos(a),sin(a)], color=red, linestyle=dash), a in angles),         seq(line([0,sin(a)],[cos(a),sin(a)], color=gray), a in {Pi/6, Pi/4, Pi/3}),         seq(line([cos(a),0],[cos(a),sin(a)], color=gray), a in {Pi/6, Pi/4, Pi/3}),         textplot([1.2,0,typeset(", ", 2*Pi)], align=right),         textplot([[1/2,0,1/2],[sqrt(2)/2,0,sqrt(2)/`2`],[sqrt(3)/2,0,sqrt(3)/`2`]], align=below),         textplot([[0,1/2,1/2],[0,sqrt(2)/2,sqrt(2)/2],[0,sqrt(3)/2,sqrt(3)/2]], align=left), axes=none, scaling=constrained, size=[800,800]);

## Use "select"...

```restart;
A := {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101};
select(x -> is(modp(x,3)=1), A);
```

`{7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97}`

By the way, if you intend A to represent the list of primes up to 101, then A should include 2 as well.  Two is the only even prime.

## Start with this and then build upon it....

As I noted earlier, it is a good idea to start with a simple case and build upon it after you gain confidence.  This toy example should be enough to get you started.

 > restart;
 > n := 2;

 > Z := Array(0..n-1,0..n-1, (i,j) -> z[i,j](t));

I have defined the coefficients arbitrarily.  Change as needed.

 > a := (i,j,k,l) -> 1*i + 2*j + 3*k + 4*l;

The (i,j) entry of the matrix on the right-hand side of the system of differential equations

 > add(add(a(i,j,k,l)*z[k,l](t), k=0..n-1), l=0..n-1): ij_entry := unapply(%, i, j);

Here is the right-hand side of the system of differential equations:

 > the_rhs := Array(0..n-1,0..n-1, ij_entry);

The system of differential equations:

 > DE := I*diff(Z, t) = the_rhs;

Initial conditions

 > IC := eval(Z, t=0) =         Array(0..n-1, 0..n-1, (i,j) -> Physics:-KroneckerDelta[i,j]);

Solve the system for the unknowns :

 > dsol := dsolve({DE, IC});

Here is the solution :

 > answer := simplify(subs(dsol, Z));

## Try this...

 > restart;
 > sigma :=  t-> (47.965 + 3.08*0.2*t^(1 - 0.7)/GAMMA(3 - 0.7))*epsilon;

Not much happens between  and :

 > plot([sigma(1), sigma(100)], epsilon=0..1);

To see a noticeable change in time, look for the time interval of zero to a million.
Here is what you get if you look at uniformly spaced time intervals:

 > tvals := seq(i, i=0..1e6, 1e5);

 > plot([seq(sigma(t), t in tvals)], epsilon=0..1);

And here is what you get with a variable time spacing:

 > tvals := seq(i^(1/(1-0.7)), i=0..100, 10);

 > plot([seq(sigma(t), t in tvals)], epsilon=0..1);

 >

## Use the AllSolutions option to int...

 > restart;
 > phi:= (mu,Q2)->sqrt(2/l)*sin(mu*Pi*(Q2+l/2)/l);

 > A := diff(phi(mu,Q2),Q2); B := diff(phi(nu,Q2),Q2);

 > -1/2/m__2*int(A*B, Q2=-l/2..l/2, AllSolutions)     assuming mu::posint, nu::posint: fh1 := combine(%);