## 6849 Reputation

16 years, 154 days

## Solution via finite differences...

This site is on the blink again and won't let me display my worksheet.   Here is the link to a worksheet that shows how to solve the system of PDEs through finite differences.

## Use this model as needed...

 > restart;
 > with(plots):
 > A1 := proc(t)         pointplot([t,t], symbol=solidcircle, symbolsize=50, color=red); end proc:
 > A2 := proc(t)         pointplot([t,t], symbol=solidcircle, symbolsize=50, color=blue); end proc:
 > display([         animate(A1, [t], t=-1..0),         animate(A2, [t], t=0..1) ], insequence);

## A bug in pdsolve()...

Your equations are uncoupled.  You may call pdsolve() to solve each of them separately.  That works on all versions of Maple that I have tried.

In principle, pdsolve() should be able to solve the two (uncoupled) equations together as a system, but that fails in Maple 2022 as we see in the attached worksheet.

Here is how things work on solving the system in various versions of Maple:

• Maple 2017 and 2018 return no solutions, and no errors;

• Maple 2019, 2020, 2021 return the correct solution;

• Maple 2022 fails on error.

 > restart;
 > kernelopts(version);

 > Physics:-Version();

Pdsolve gets the correct solution here:

 > pde1 := diff(u(x,t),t) = 0; ic1 := u(x,0) = f(x);

 > pdsolve({pde1,ic1});

Pdsolve gets also the correct solution here:

 > pde2 := diff(v(x,t),t) = 0; ic2 := v(x,0) = g(x);

 > pdsolve({pde2,ic2});

BUG: Pdsolve gets confused when solving the two (uncoupled) equations together:

 > pdsolve({pde1, pde2, ic1, ic2});

Error, (in pdsolve) invalid input: indets expects 1 or 2 arguments, but received 3

 >

## Getting started...

I am not going to do the entire homework for you but here is something to get you started.

 > plot(8+8*sin(theta), theta=0..2*Pi, coords=polar, scaling=constrained);
 > plot3d(25 - x^2, x=0..5, y=1..8, view=[0..5,0..8,default]);
 > plots:-display(     plot3d([10*cos(t)*cos(s), 10*cos(t)*sin(s), 10*sin(t)], t=-arccos(6/10)..arccos(6/10), s=-Pi..Pi),     plot3d([6*cos(t), 6*sin(t), s], t=-Pi..Pi, s=-8..8), color=["Green","Red"], style=surface, scaling=constrained, lightmodel=light4);

## A suggestion for an alternative...

I have not examined the source of the error in your worksheet, but I would suggest using the facilities provided in the Physics[Vectors] package for encoding the Rodrigues formula rather than homogeneous transformation matrices, because the math works more naturally this way.  This worksheet shows how.

 > restart;
 > with(Physics[Vectors]):

Rodrigues'  formula:
Rotates the vector about the unit vecor  by the angle :

 > v__rot_ := v_*cos(alpha) + (k_ &x v_)*sin(alpha) + (k_ . v_)*(1 - cos(alpha))*k_;

Make that into a procedure: Usage:

 > R := unapply(v__rot_, v_, k_, alpha);

Example 1

Pick a unit vector:

 > 1*_i + 2*_j + 3*_k: k_ := % / sqrt(% . %);

Pick an arbitrary vector:

 > v_ := _i - 3*_j + 2*_k;

Pick a rotation angle:

 > alpha := Pi/6;

Find the rotated vector:

 > R(v_, k_, alpha): collect(%, [_i, _j, _k]);

 >

Example 2

Composite rotation

A unit vector:

 > k__1_ := k_;

Another unit vector

 > k__2_ := subs(_j=-_j, k__1_);

 > R(R(v_, k__1_, Pi/6), k__2_, Pi/4): collect(%, [_i, _j, _k], simplify);

 >
 >

## Solution...

 > restart;
 > kernelopts(version);

 > pd1 := diff(u(x,t),t) = diff(u(x,t),x,x) + (1-alpha)*u(x,t);

 > pd2 := diff(v(x,t),t) = mu*diff(v(x,t),x,x) + beta*v(x,t) + alpha*u(x,t);

 > ics := u(x,0) = Dirac(x), v(x,0) = 0;

Maple is unable solve the system sympolically on its own:

 > pdsolve({pd1, pd2, ics}) assuming t > 0, mu > 0;

but it can, with a little help.

The first PDE has only the unknown , so we solve it for :

 > sol1 := pdsolve({pd1,ics[1]}) assuming t > 0;

Then, substitute the result in the second PDE, and solve that for :

 > subs(sol1, pd2); sol2 := pdsolve({%, ics[2]}) assuming t > 0, mu > 0;

## Braces, brackets, and parentheses...

In several places you have used square brackets [...] and curly braces {...} for grouping your mathematical terms.  In Maple, the only delimiters that you can use for that purpose are the parentheses (...).

Adjust your code as necessary and write back if it still does not work.

## Simplify...

delvin, you should know that your equations are not comprehensible to people who are unfamiliar with numerals written in the Farsi/Persian script.  If you desire to reach the widest available help, make an effort to write your mathematics in symbols which most of the readers of this forum can understand.

 > restart;

You wish to perform your hand calculation in Maple.  Here is the original expression:

 > c := (-exp(4*d)*beta^2 + 4*alpha*exp(4*d) + beta^2 - 4*alpha)/(8*exp(2*d));

and here is the simplified form obtained by Maple:

 > convert(simplify(c), trigh);

That agrees with your hand calculation.

## This ought to do it...

 > restart;
 > g := x -> piecewise(x - 1/3 = 0, 1, 0);

 > g(3+1/3);

 > g(0+1/3);

## Solution...

 > restart;

Take  and renane  to :

 > 1/2*y^2 + (1-cos(phi)) = h+2: solve(%, y): f := unapply(%[1], phi, h);

 > seq(f(phi,h), h in [-1,0,1]), seq(-f(phi,h), h in [-1,0,1]); plot([%], phi=-Pi..3*Pi, size=[800,500]);

## Infinitely many solutions...

That initial value problem has infinitely many solutions.  Have a look at this worksheet.

 > restart;
 > ode := diff(v(t),t) = -2*v(t)^(2/3);

Find the general solution

 > dsolve(ode); dsol := isolate(%, v(t));

Let :

 > eval(dsol, c__1 = 0);

That's a solution of the ODE with the initial condition .

That, however, is not the only solution.  The following is also

a solution of the initial value problem for any choice of

 > sol := piecewise(t < a, 0, -(8*(t-a)^3)/27);

Here is what the solution looks like when :

 > eval(sol, a=2); plot(%, t=0..4, color=red, thickness=5);

## Assumed variables...

The trailing tilde is not a part of  Omega; it is there to remind you that you have made assumptions on it.

If you don't want to see those reminders, disable them with inserting interface(showassumed=0); near the top of your worksheet.

## Convert to D...

mtaylor does not like expressions with "diff" in them, but it works well when converted to "D".  So all you need is:

```convert(lhs(Eq1), D):
mtaylor(%, [x,t], 2);```

## Geodesics between two points on a torus...

Geodesics connecting two points on a torus

 > restart;
 > with(plots):
 > with(VariationalCalculus):

The parametric equation of a torus with major and minor radii of  and .

 > < (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >; x := unapply(%, u, v):

Calculate the coefficients , , and  of the first fundamental form:

 > xu := diff(x(u,v),u): xv := diff(x(u,v),v): E := unapply(simplify(xu^+ . xu), u, v); F := unapply(simplify(xu^+ . xv), u, v); G := unapply(simplify(xv^+ . xv), u, v);

Look for geodesics of the form .  These account for all geodesics with the exception

of those corresponding to constant.  The latter are plain circles so we don't need Maple

to calculate them.

From differential geometry, the length of  any curve of the form

is given by  where:

 > Upsilon := sqrt(E(u, v(u)) + 2*F(u, v(u))*diff(v(u),u) + G(u,v(u))*(diff(v(u),u))^2);

The minimum length is obtained with the help from the Euler-Lagrange equation

which leads to a second order nonlinear differential equation in :

 > EulerLagrange(Upsilon, u, v(u)): remove(type, %, equation)[]: isolate(%, diff(v(u), u, u)): de := simplify(%);

Pick numbers for the torus's radii:

 > a, b := 5, 2;

 > torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi,         scaling=constrained, color="khaki", style=surface);

Example 1

Calculate a geodesic from  to ,  :

 > bc1 := v(0)=0,  v(Pi/2)=Pi/2;

 > dsol1 := dsolve({de, bc1}, numeric, output=operator);

 > display(         torus,         tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),         pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol1), symbol=solidsphere, symbolsize=20, color="Orange"), lightmodel=light4, thickness=5, orientation=[-5,75,0], size=[700,400], style=surface, axes=none);

 >

Example 2

 > bc2 := v(0)=0,  v(Pi/2)=Pi;

 > dsol2 := dsolve({de, bc2}, numeric, output=operator);

 > display(         torus,         tubeplot([seq(eval(x(u, v(u)), dsol2))], u=0..Pi/2, color=red, radius=0.1),         pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol2), symbol=solidsphere, symbolsize=20, color="Orange"), lightmodel=light4, thickness=5, orientation=[-65,65,0], size=[700,400], style=surface, axes=none);

 >
 >

Example 3

The end points in this example are identical to those of Example 1 but the new geodesic,
drawn in cyan, is different.  In fact, there are infinitely many  geodesics that connect

those two points.

 > bc3 := v(0)=0,  v(Pi/2)=5*Pi/2;

 > dsol3 := dsolve({de, bc3}, numeric, output=operator);

 > display(         torus,         tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),         tubeplot([seq(eval(x(u, v(u)), dsol3))], u=0..Pi/2, color=cyan, radius=0.1),         pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol3), symbol=solidsphere, symbolsize=20, color="Orange"), lightmodel=light4, thickness=5, orientation=[-65,45,0], size=[700,400], style=surface, axes=none, viewpoint=circleright);

 >
 >
 >

## Finding geodesics...

You will find the geodesic equations in the wikipedia web page that you have referenced to.   Admittedly it is rather cryptic, involving Christoffel symbols in an implicit summation.  Once we unpack the notation, the geodesic equations for a surface turn out to be a pair of coupled, second order, nonlinear ODEs.  These equations are usually derived in an introductory course on differential geometry.  To get to those equations, one needs to learn about a surface's first fundamental form and how the Christoffel symbols are expressed in terms of the coefficients of the first fundamental form.

The Christoffel() proc in the attached worksheet calculates the coefficients of the first fundamental form of any surface and then calculates and returns the corresponding Christoffel symbols.  To understand the details, you need to know some differential geometry.

A second proc, Geodesic_DEs(), applies the equation in the Wikipedia page to calculate and return the surface's geodesic equations.  Once you have those equations, you can apply Maple's dsolve() to calculate and plot the geodesics.  The worksheet includes an illustrative example.

Maple's Differential Geometry package has a GeodesicEquations() function.  Odds are that what I have calculated here can be done more quickly by applying that function but don't know how.  Perhaps someone who knows may show us.

Answers to the several questions that you asked would involve some elementary differential geometry, so I cannot summarize them here in a few lines.  I can, however, give short answers to two of those questions:

1. What does "to walk forward" mean?  When you move on a surface, at each location you have a velocity vector, an acceleration vector, and a normal vector (i.e., vector perpendicular to the surface).  You are "moving forward", or equivalently, "traveling on a geodesic", provided that those three vectors are co-planar.
2. Geodesic is the shortest distance between two points in space, but no such end points are shown in the website's animation.  A geodesic typically extends indefinitely in both directions.  There are no end points.  A curve is geodesic if for any pair of points A, B on the curve, the part of the curve that lies between A and B is the shortest path among all curves that lie on the surface and connect A to B.
 > restart;
 > with(plots):

Receives a parametrization  of a surface and returns the corresponding

Christoffel symbols , each of which is a procedure.

 > Christoffel := proc(x::procedure)         local u, v, xu, xv, E, F, G, Eu, Ev, Fu, Fv, Gu, Gv, den, C;         xu := diff(x(u,v),u);         xv := diff(x(u,v),v);         # E, F, G are the coefficients of the surface's first fundamental form         E := xu^+ . xu;         F := xu^+ . xv;         G := xv^+ . xv;         Eu := diff(E, u);         Ev := diff(E, v);         Fu := diff(F, u);         Fv := diff(F, v);         Gu := diff(G, u);         Gv := diff(G, v);         den := 2*(E*G - F^2);         C[1][1,1] := unapply(simplify((G*Eu - 2*F*Fu + F*Ev)/den), u, v);         C[1][1,2] := unapply(simplify((G*Ev - F*Gu)/den), u, v);         C[1][2,2] := unapply(simplify(-(F*Gv - 2*G*Fv + G*Gu)/den), u, v);         C[2][1,1] := unapply(simplify(-(E*Ev - 2*E*Fu + F*Eu)/den), u, v);         C[2][1,2] := unapply(simplify((E*Gu - F*Ev)/den), u, v);         C[2][2,2] := unapply(simplify((E*Gv - 2*F*Fv + F*Gu)/den), u, v);         C[1][2,1] := C[1][1,2];         C[2][2,1] := C[2][1,2];         return C; end proc:

Example: A sphere

 > x := (u,v) -> < cos(v)*cos(u), cos(v)*sin(u), sin(v) >;

Here are the 8 Christoffel symbols of our parametrization of the sphere:

 > Gamma := Christoffel(x):
 > Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

 > Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

 > x := (u,v) -> < (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >;

 > Gamma := Christoffel(x):
 > Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

 > Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

The geodesic equations

 > Geodesic_DEs := proc(x::procedure)         local Gamma, de, i, j, k, X;         Gamma := Christoffel(x);         for i from 1 to 2 do            de[i] := diff(X[i](t), t, t) + add(add(                         Gamma[i][j,k](X[1](t),X[2](t))*diff(X[j](t),t)*diff(X[k](t),t),                         j=1..2), k=1..2) = 0;         end do:         eval([de[1],de[2]], {X[1]=u, X[2]=v})[];   # change notation from X to u, v         return (%); end proc:
 >

Geodesics on a torus

Here we calculate and plot geodesics on a torus.  You may change the torus

to any other parametrized surface.  Everthing else remains the same.

 > < (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >; eval(%, {a=5, b=2}); x := unapply(%, u, v):

 > the_torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi, scaling=constrained, style=wireframe);

The differential equations of geodesics on torus

 > DEs := Geodesic_DEs(x);

Pick any desired starting point and direction on the torus.  The angle  determines
the geodesic's initial direction.

 > unassign('theta'); ICs := u(0) = 0, v(0) = 0, D(u)(0) = cos(theta), D(v)(0) = sin(theta);

Pick a  and solve the geodesic DEs

 > theta:= 0.95*Pi/2: dsol := dsolve({DEs, ICs}, numeric, output=operator);