Rouben Rostamian

Symbolic solution...

Answered: Rouben Rostamian 8636

September 03 2023

4 0

Maple can find a symbolic solution to your boundary value problem in terms of special functions. There is no need to resort to a numerical approximation.

PS: Note added later:

It appears that you have imposed the odd-looking boundary condition y(1e-14)=1 solely to avoid division by zero in your numerical calculation. The symbolic solution shows that the limit of the solution as x approaches zero from the right is 1, so you might as well change your boundary condition to y(0)=1. To verify that, set eps=0 in the worksheet that I posted earlier.

>

restart;

>	Digits := 30;

>	eps := 1.0e-14;

>	de := sqrt(x)diff(y(x),x,x) - 3/2y(x) - 1/2 = 0;

x^(1/2)*(diff(diff(y(x), x), x))-(3/2)*y(x)-1/2 = 0

>	dsol := dsolve({de, y(eps)=1, y(10)=0});

>	indets(dsol, name);

>	evalf(subs(x=eps, dsol));

>	evalf(subs(x=10, dsol));

>	plot(rhs(dsol), x=0..10);

>

Download zz.mw

A numerical solution...

Answered: Rouben Rostamian 8636

August 26 2023

1 7

If you do:

fsolve({Eq1=p_o-p_i,Eq2=F_a}, {A=0.167629..0.167634, B=1.000640 .. 1.000650});

You will get:

{A = 0.1676309260049029044111806984482988970179873516977323,
 B = 1.000644958082968954903634017487303572859626170659325}

I obtained the rough range of the roots by trial-and-error. There may be other roots.

Solution...

Answered: Rouben Rostamian 8636

August 24 2023

2 2

Calculating the tangent plane to a sphere is quite trivial. Why bother with the geometry package for it?

Just write a proc to receive a point on the sphere, calculate the tangent plane, and format the result according to your specification. Then map the proc over your list of points.

>

restart;

>	S := (x-1)^2 + (y-3)^2 + (z-5)^2 = 13^2;

>

pts := [[-12, 3, 5], [-11, -2, 5], [-11, -1, 2], [-11, -1, 8], [-11, 0,
        1], [-11, 0, 9], [-11, 3, 0], [-11, 3, 10], [-11, 6, 1], [-11, 6,
        9], [-11, 7, 2], [-11, 7, 8], [-11, 8, 5], [-4, -9, 5], [-4,
        3, -7], [-4, 3, 17], [-4, 15, 5], [-3, -9, 2], [-3, -9, 8], [-3,
        0, -7], [-3, 0, 17], [-3, 6, -7], [-3, 6, 17], [-3, 15, 2], [-3, 15,
        8], [-2, -9, 1], [-2, -9, 9], [-2, -1, -7], [-2, -1, 17], [-2,
        7, -7], [-2, 7, 17], [-2, 15, 1], [-2, 15, 9], [1, -10, 5], [1, -9,
        0], [1, -9, 10], [1, -2, -7], [1, -2, 17], [1, 3, -8], [1, 3,
        18], [1, 8, -7], [1, 8, 17], [1, 15, 0], [1, 15, 10], [1, 16,
        5], [4, -9, 1], [4, -9, 9], [4, -1, -7], [4, -1, 17], [4,
        7, -7], [4, 7, 17], [4, 15, 1], [4, 15, 9], [5, -9, 2], [5, -9,
        8], [5, 0, -7], [5, 0, 17], [5, 6, -7], [5, 6, 17], [5, 15, 2], [5,
        15, 8], [6, -9, 5], [6, 3, -7], [6, 3, 17], [6, 15, 5], [13, -2,
        5], [13, -1, 2], [13, -1, 8], [13, 0, 1], [13, 0, 9], [13, 3,
        0], [13, 3, 10], [13, 6, 1], [13, 6, 9], [13, 7, 2], [13, 7,
        8], [13, 8, 5], [14, 3, 5]]:

>

tangent_plane := proc(P)
        local C, N, ans;
        C := [1,3,5];      # the center of the sphere
        N := P - C;        # the normal vector at P
        (x - P[1])*N[1] + (y - P[2])*N[2] + (z - P[3])*N[3]; # <x,y,z> . N
        # divide the coefficients by their greatest common divisor
        igcd(coeffs(%));
        # we want the sign of the leading coefficient to be positive
        ans := %%/%;
        if coeff(%, x) > 0 then
                return [P,ans = 0];
        elif coeff(%, x) < 0 then
                return [P,-ans = 0];
        elif coeff(%, y) > 0 then
                return [P,ans = 0];
        elif coeff(%, y) < 0 then
                return [P,-ans = 0];
        elif coeff(%, z) >= 0 then
                return [P,ans = 0];
        else
                return [P,-ans = 0];
        end if;
end proc:

>	map(tangent_plane, pts);

Download mw.mw

Animate pursuit...

Answered: Rouben Rostamian 8636

August 18 2023

2 1

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.

Download pursuit.mw

No solution?...

Answered: Rouben Rostamian 8636

August 14 2023

1 2

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. ...

Answered: Rouben Rostamian 8636

August 09 2023

0 2

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;

Download New-eval-2.mw

Built-in...

Answered: Rouben Rostamian 8636

August 01 2023

1 3

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

Download mw.mw

A solution...

Answered: Rouben Rostamian 8636

July 31 2023

0 1

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 := <p, acos(s), asin(s) >;

Vector(3, {(1) = p, (2) = a*cos(s), (3) = a*sin(s)})

Cylinder of radius centered on the axis; parameters and

>	C2 := <ccos(t), csin(t), q>;

Vector(3, {(1) = c*cos(t), (2) = c*sin(t), (3) = q})

Intersection curve

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

Vector(3, {(1) = p = c*cos(t), (2) = a*cos(s) = c*sin(t), (3) = a*sin(s) = q})

${p = c*cos(t), q = a*sqrt(-(c^2*sin(t)^2-a^2)/a^2), s = arccos(c*sin(t)/a)}$

Vector[column](%id = 36893628579494067244)

>

%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, %));

>

Download intersection-of-cylinders.mw

Elliptic section...

Answered: Rouben Rostamian 8636

July 18 2023

3 0

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>;

Vector(3, {(1) = s*tan(alpha)*cos(t), (2) = s*tan(alpha)*sin(t), (3) = s})

>	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);

Vector(3, {(1) = sin(beta), (2) = 0, (3) = cos(beta)})

Vector(3, {(1) = 0, (2) = 1, (3) = 0})

Vector[column](%id = 36893628788357808124)

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

>	P := Au + Bv + d*<0,0,1>;

Vector(3, {(1) = -v*cos(beta), (2) = u, (3) = v*sin(beta)+d})

>	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);

Vector(3, {(1) = d*cos(beta)*tan(alpha)*cos(t)/(sin(beta)*tan(alpha)*cos(t)+cos(beta)), (2) = d*tan(alpha)*sin(t)*cos(beta)/(sin(beta)*tan(alpha)*cos(t)+cos(beta)), (3) = d*cos(beta)/(sin(beta)*tan(alpha)*cos(t)+cos(beta))})

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

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

Download elliptic-section.mw

Try this simpler one...

Answered: Rouben Rostamian 8636

July 15 2023

6 4

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);

Matrix(4, 4, {(1, 1) = 50, (1, 2) = -50, (1, 3) = -38, (1, 4) = -98, (2, 1) = 10, (2, 2) = -22, (2, 3) = -18, (2, 4) = -77, (3, 1) = -16, (3, 2) = 45, (3, 3) = 87, (3, 4) = 57, (4, 1) = -9, (4, 2) = -81, (4, 3) = 33, (4, 4) = 27})

(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);

Matrix(4, 4, {(1, 1) = 50., (1, 2) = 1445.5341186930428984, (1, 3) = -38., (1, 4) = -98., (2, 1) = 10., (2, 2) = -22., (2, 3) = -18., (2, 4) = -77., (3, 1) = -16., (3, 2) = 45., (3, 3) = 87., (3, 4) = 57., (4, 1) = -9., (4, 2) = -81., (4, 3) = 33., (4, 4) = 27.})

(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]^+;

Vector(4, {(1) = .986301870755100907032581022014, (2) = -0.236750758506050268946221115509e-1, (3) = .124321775276849024311433861988, (4) = .105793226250403080836382592166})

(9)

>	M . null__vec;

Vector(4, {(1) = 0.2064116200e-20, (2) = 0.2275617400e-21, (3) = -0.1073762554e-19, (4) = 0.3080927157e-19})

(10)

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

that is in the null space.

Download singular_values.mw

Corrections...

Answered: Rouben Rostamian 8636

July 12 2023

5 6

>

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)};

{RR(2) = 2, diff(RR(r), r) = 1/(1-1/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
dq/dr = 1/(1-1/r) 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.

Download mw.mw

The proper demo...

Answered: Rouben Rostamian 8636

June 30 2023

3 2

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...

Answered: Rouben Rostamian 8636

June 23 2023

2 2

>	restart; q1:=1-1/(w*2-3sigmak2)-deltabmu/((w-ku0b)2-3mudeltabsigmabk2)+A/k*2=0;

1-1/(-3*k^2*sigma+w^2)-deltab*mu/((-k*u0b+w)^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);

1-1/(Q^2*k^2-3*k^2*sigma)-deltab*mu/((Q*k-k*u0b)^2-3*mu*deltab*sigmab*k^2)+A/k^2 = 0

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);

1-1/(Q^2*k^2-3*k^2)-mu/((Q*k-k)^2-3*mu*k^2)+1/k^2 = 0

>	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.

Download mw.mw

A clean way...

Answered: Rouben Rostamian 8636

June 21 2023

1 1

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(iPi/4, i=0..7), seq(iPi/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]);

Download mw.mw

Use "select"...

Answered: Rouben Rostamian 8636

June 20 2023

1 1

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);

The answer is

{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.

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

8636 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Symbolic solution...

A numerical solution...

Solution...

Animate pursuit...

No solution?...

See if this does what you want. ...

Built-in...

A solution...

Elliptic section...

Try this simpler one...

Corrections...

The proper demo...

There is no neat solution...

A clean way...

Use "select"...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again