## 7776 Reputation

17 years, 361 days

## Use quotes...

It works if you place h(p) in quotes.  I don't know why.

`plot('h(p)', p=-1..1);`

Doing

`showstat(MapleTA:-Builtin:-inverf)`

we see that inverf() is defined as:

`inverf := proc(z::Range(0,1)) ...`

and therefore it accepts only numerical arguments in the range (0,1).  In particular, inverf(p) is not accepted if p is an unassigned symbol.  That's bad.  A better-designed inverf() would have returned inverf(p) (unevaluated) for such a p, and that would have made the quotation marks in the plot() command unnecessary.

## Procedure for calculating the Christoffe...

In this worksheet I show how to calculate the Christoffel symbols for a two-dimensional manifold, that is, a surface embedded in 3D.

Christoffel symbols of a parametrized surface

 > restart;

Receives a metric tensor as a  matrix in the variables  and .

Returns the Christoffel tensor  which is represented as

Gamma[k][i,j] in Maple.

The names  are not hard-coded.  They can be anything.

 > Christoffel := proc(g::Matrix(2,2), u::assignable, v::assignable)         local U, G, C, i, j, k, l;         U[1], U[2] := u, v;         G := g^(-1);         for i from 1 to 2 do                 for j from 1 to 2 do                         for l from 1 to 2 do                                 C[l][i,j] := 1/2*add(G[k,l]*(diff(g[i,k], U[j]) - diff(g[i,j], U[k]) + diff(g[k,j], U[i])), k=1..2);                         end do;                 end do;         end do;         return C; end proc:
 >

Examples of use

 > S := ;

 > plot3d(S, x=-1..1, y=-1..1);

 > X[1] := diff(S, x): X[2] := diff(S, y): g := simplify(Matrix(2,2, (i,j) -> X[i]^+ . X[j]));

 > Gamma := Christoffel(g, x, y):

Here are the 8 components of :

 > simplify([Gamma[1][1,1], Gamma[1][1,2], Gamma[1][2,1], Gamma[1][2,2]]); simplify([Gamma[2][1,1], Gamma[2][1,2], Gamma[2][2,1], Gamma[2][2,2]]);

A helicoid

 > S := < s*cos(t), s*sin(t), t/4>;

 > plot3d(S, t=0..4*Pi, s=0..1, grid=[80,8]);

 > X[1] := diff(S,t): X[2] := diff(S,s): g := simplify(Matrix(2,2, (i,j) -> X[i]^+ . X[j]));

Here we write  for the Christoffel symbols instead of the usual  to show that

the name  is not special:

 > K := Christoffel(g, t, s):
 > K[1][1,1], K[1][1,2], K[1][2,1], K[1][2,2]; K[2][1,1], K[2][1,2], K[2][2,1], K[2][2,2];

A torus

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

 > plot3d(subs(a=5, b=2, S), u=-Pi..Pi, v=-Pi..Pi, scaling=constrained);

 > X[1] := diff(S,u): X[2] := diff(S,v): g := simplify(Matrix(2,2, (i,j) -> X[i]^+ . X[j]));

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

A generic metric tensor

 > g := < a(u,v), b(u,v); b(u,v), c(u,v) >;

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

A generic parametric surface

 > S := < alpha(s,t), beta(s,t), gamma(s,t) >;

 > X[1] := diff(S,s): X[2] := diff(S,t): g := simplify(Matrix(2,2, (i,j) -> X[i]^+ . X[j]));

 > Gamma := Christoffel(g, s, t):
 > simplify([Gamma[1][1,1], Gamma[1][1,2], Gamma[1][2,1], Gamma[1][2,2]]); simplify([Gamma[2][1,1], Gamma[2][1,2], Gamma[2][2,1], Gamma[2][2,2]]);

 >

## plotsetup...

If you want your plot to appear in a separate window, do one or the other of the following options;

```# plot in a separate window
plotsetup(window);
plot(x^2, x=-1..1);

# plot as a maplet
plotsetup(maplet);
plot(x^2, x=-1..1);

# restore the default (inline) plotting mode
plotsetup(default);
plot(x^2, x=-1..1);
```

## See if this helps...

 > restart;
 > n := 5;

 > LinearAlgebra:-BandMatrix([-1, 2, -1], 1, n);

## Number of roots = 116 is correct...

 > restart;
 > with(plots):
 > setoptions3d(style=surface);
 > f1 := x3^2-1/10*x1^4-5/100*x2^4+1; f2 := x1^3+x2^3+5/100*x3^3-1; f3 := -2*cos(3*x1)+2*cos(3*x2)-2*cos(3*x3)+1;

Relevant ranges of x1, x2, x3, determined with plotting:

 > a, b, c := 10, 14, 35;

 > implicitplot3d(f1, x1=-a..a, x2=-b..b, x3=-c..c); p1 := %:

 > implicitplot3d(f2, x1=-a..a, x2=-b..b, x3=-c..c); p2 := %:

 > implicitplot3d(f3, x1=-a..a, x2=-b..b, x3=-c..c); p3 := %:

Want to determine the equation of the intersection curve of the surfaces f1 and f2:

 > implicitplot3d([f1,f2], x1=-a..a, x2=-b..b, x3=-c..c,         color=[red,blue], lightmodel=light4, orientation=[40,80,0]);

To find the intersection, parametrize the intersection curve as  and find a differential
equation for it.   Append the equation  to make the parametrization into

a unit-speed curve:

 > subs(x1=x(t), x2=y(t), x3=z(t), [f1,f2]): diff(%, t): EQs := %[], diff(x(t),t)^2 + diff(y(t),t)^2 + diff(z(t),t)^2 - 1 ;

Isolate the derivatives in the equations above.  We get two solutions, corresponding to traversing

the curve in one or the other directions.  Without loss of generality pick one of the two options.

 > solve({EQs}, {diff(x(t),t), diff(y(t),t), diff(z(t),t)}): allvalues(%): DEs := %[1];

Now we need an initial condition for our set of DEs.   Any point on the intersection curve will do.
So we set x3=20 arbitrarily and solve for x1 and x2;

 > subs(x3=20, [f1, f2]): fsolve(%); ICs := subs(%, {x(0)=x1, y(0)=x2, z(0)=20});

We are ready to solve the system of DEs:

 > dsol := dsolve(DEs union ICs, numeric, output=operator);

Here is what the intersection curve looks like.  The range (that is, the curve's period)

is determined by trial-and-error so that to have a closed curve:

 > eval([x(t),y(t),z(t)], dsol): plots:-tubeplot(%, t=0..180, color=yellow, radius=0.5,         numpoints=200, orientation=[25,75,0], labels=[x,y,z]); curve := %:

Let's verify that we have indeed found the correct intersection;

 > display(p1,p2,curve, color=[red,blue,yellow], orientation=[25,75,0]);

Finally, we want to find the intersection points of the curve we have just

found with the surface f3=0.  So we evaluate f3 along the curve over one

period and count the number of zeros:

 > subs(x1=x(t), x2=y(t), x3=z(t), f3): eval(%, dsol): plot(%, t=0..180); zeros := fsolve(%%, t=0..180, maxsols=300);

How many zeros?

 > nops([zeros]);

 >

## An elementary solution...

Let's write a for the length of AD.  The fact that we are given a=4 is not essential.  I will write q for the angle ACD, and x for the length of AB.

Apply the law of sines in the triangle ACD:
eq1 := sin(80) / c = sin(q) / a     # forgive me for using degrees for the argument of sin()

Apply the law of sines in the triangle ACB:
eq2 := sin(80+20) / (c + 4/3*c) = sin(q) / x.

Dividing the respective sides of eq1 and eq2, we obtain:
7/3*sin(80) / sin(100) = x / a.

We note, however, that sin(80) = sin(90-10) = sin(90+10) = sin(100), and therefore the previous equation reduces to

x / a = 7 / 3.

In particular, if a = 4, we get x = 28 / 3.

I have used Maple in nothing but Linux for longer than I can remember and haven't had problems such as those that you have described.  Perhaps your Linux installation is broken?

That said, it is possible to run out of java heap memory if you are doing heavy-duty graphics.  In that case consider specifying a larger java heap on the command line, as in

```maple -x -j 65536 &
```

The default java heap is 512.

PS: My current Linux is Ubuntu 22.04.

I see several issues with your worksheet.

• You have three ordinary (not partial) differential equations in e1, e2, e3.  The command to solve ODEs is dsolve(), not pdsolve().
• It looks like the only unknown in your equations is z2(t).  Therefore each of e1, e2, e3 may be solved independently.  There is no guarantee that the z2(t) produced by one of these equations would agree with the z2(t) produced by the other.  Why are you solving the equations as a system?
• Each of your e1, e2, e3 involves a coefficient _C1.  Naming a symbol with a leading underscore is a bad idea.  Such symbols are reserved for use by Maple.  Your _C1 can potentially conflict with Maple's operation.

## Solution...

There may or may not be a solution, depending on some assumptions which you have not stated.

For instance, if you assume that x1_prime, x2_prime, and u are constants, then a solution exists and can be obtained through

```pde1 := x1_prime = diff(H(x1,x2),x1) + u*x1 + x2;
pde2 := x2_prime = x1 + diff(H(x1,x2),x2);
pdsolve({pde1,pde2});```

## If you are comfortable with C programmin...

If you are comfortable with C programming, you should be able to write your own implementation of KroneckerProduct() in C and compile it with the rest of your program.  When you call f() as exported by Maple, when it hits the call to KroneckerProduct(), it will execute your own version of that function.

The effective writing of such a KroneckerProduct() function will depend on your C programming skills, and some details that are missing in what you have written.  Here are some thoughts:

1. In your example, A, B, and their Kronecker product are single-column matrices.  If single-columns matrices are all you need, then the C program becomes significantly simpler because you may use one-dimensional arrays to represent them.
2. In your example, A and B have 3 rows each.  If 3-row matrices are all you need, then the C program becomes even simpler.
3. How familiar are you with the various C standards?  The C99 standard introduced variable length arrays, which can help with your program.  The previous standard, C89, recognizes only arrays whose sizes are known at the compile time.  Variable length arrays may be simulated in C89 through dynamic memory allocation with malloc() and related functions.  You need to decide which method you want to use.

## Needs better care...

• Look at your s[2].  On the left, you have E(t).  On the right you have E.  That should be E(t).  Also on the right you have lambda[2] and lambda[3].  Those should be lambda[2](t) and lambda[3](t).  Similar errors occurs in several other places.  Need to go over your input very carefully and fix those things.
• In s[5] you have lambda[1]*t.  That should be lambda[1](t).
• You have set tf := T, but you didn't say what T is.  Need to specify a numerical value for T.
• Don't apply Maple commands blindly.  You have defined fcns[1] := ... and ended that line with a colon.  Replace the colon with a semicolon to verify that Maple understands you correctly.  Similarly, you call dsolve with the arugment
`{seq(s[i], i = 1 .. 8), seq(lambda[i](tf) = 0, i = 1 .. 4), E(0) = E0, In(0) = In0, R(0) = R0, S(0) = S0}`

Before calling dsolve, enter that argument on a line of its own and examine what Maple prints very carefully to verify that what you are later passing to dsolve is in agreement with what you have in mind.

## Maybe this?...

I am not familiar with Mathematica, so I am just guessing from your example what the Cases command is expected to do.  It appears that it picks all instances of sin(...) in a given expression.  Here is how that is done in Maple:

expr := [sin(x)/(sin(2-x)+1)-12, sin(x/2)^2, cos(x)];
indets(expr, 'specfunc(sin)');

That yields {sin(x), sin(-2 + x), sin(x/2)}.

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

 >