## 7776 Reputation

17 years, 361 days

## Remove the quotes...

The construction `A=70` (with the back-quotes) builds a name consisting of those four characters, just like AB70 is a name consisting of four characters.  The equal sign has no special meaning in that name.  If instead you do

`lst := [A=70, B=17, C=27];`

(without the quotes,) then

`tt := convert(lst, table);`

will do what you want.

## Ill-posed problem...

This problem is not well-posed.  See the discussion here.

## That may not be a meaningful operation...

 What is the meaning of ? The Dirac delta is defined through the property     for all continuous functions .   Let's see if we can make sense of something simpler than your exmaple, say . Make a change of variables  and calculate . To continue, we need to evaluate the expression within the square brackets at  which may not be feasible for all  and therefore  makes no sense.

## Will this do?...

```restart;
with(plots):
with(geom3d):
ball := draw(TruncatedIcosahedron(``)(color="Green"));
```

## Series calculation...

This worksheet shows how to calculate the first five terms of the desired series.  I have coded it as 5 similar blocks.  They may be combined into single for-loop to compute any number of terms.

The coefficient of u^3 in what you have shown seems to be in error.  The (a-1)(a+3) should be (a-1)(a+2)..

We have

This defines  as a function of  but we want  as a function of ,

as in  so we calculate the Maclaurin series of .
Observe that  implies that , and therefore

.  That gives us the first term of the series.

We find the remaining terms by repeated differentiations.

 > restart;
 > eq[0] := z = int(sinh(a*sin(x(z)*y)), y=0..1);

 > rels[0] := x(0)=0;

 > eq[1] := diff(eq[0], z); convert(%, D); eval(%, z=0); eval(%, {rels[0]}); isolate(%, D(x)(0)); rels[1] := rels[0], simplify(%);

 > eq[2] := diff(eq[1], z): convert(%, D); eval(%, z=0); eval(%, {rels[1]}); isolate(%, (D@@2)(x)(0)); rels[2] := rels[1], simplify(%);

 > eq[3] := diff(eq[2], z): convert(%, D); eval(%, z=0); eval(%, {rels[2]}); isolate(%, (D@@3)(x)(0)); rels[3] := rels[2], simplify(%);

 > eq[4] := diff(eq[3], z): convert(%, D); eval(%, z=0); eval(%, {rels[3]}); isolate(%, (D@@4)(x)(0)); rels[4] := rels[3], simplify(%);

 > eq[5] := diff(eq[4], z): convert(%, D); eval(%, z=0); eval(%, {rels[4]}); isolate(%, (D@@5)(x)(0)); rels[5] := rels[4], simplify(%);

Here is the series:

 > series(x(z), z); eval(%, {rels[5]});

 >

## Errors...

You have four PDEs in the five unknowns

`C(X, R, t), R(X, R, t), T(X, R, t), U(X, R, t), V(X, R, t)`

Additionally, U(X,R,t) appears just as U in the fourth PDE.  There may be other errors as well.

Suggestions:

1. Post an actual worksheet rather than copy/paste.  That will help people that are going to help you.
2. Before you do a composite computation in a for-loop, try just one case in order to detect where the error occurs.  You may add the for-loop afterward once you see that the basic code works.

## Is this good enough?...

```restart;
ode1 := diff(x(t),t)=2*x(t)*y(t);
ode2 := diff(y(t),t)=1-x(t)^2-y(t)^2;
DEtools:-DEplot([ode1,ode2],[x(t),y(t)],
t=0..1,x = -4..4, y = -4..4,
arrows=curve, linecolor=red, arrowsize=1.5,
axes=boxed, color='magnitude[legacy]');
```

## Heap size...

The default java heap size is 512.  You may change it through the commandline flag −j, as in

`maple -j 65536 -x filename.mw`

## Use "Explore"...

Use Explore.  Move the sliders to vary M and k.

 > restart;
 > addcoords(z_cylindrical, [z, r, t], [r*cos(t), r*sin(t),z]):
 > F := -2*M^2*sin(theta)^2 - 2*k*r^2

 > Explore(plot3d([F,0], r=0..1, theta=0..2*Pi, coords=z_cylindrical), M=-1.0..1.0, k=-1.0..1.0);

## Maybe this?...

It's hard to understand what you are asking.  Maybe this is what you are looking for?

 > restart;
 > with(plots):

Two space curves, parametrized as :

 > y := sqrt(1 - x^2); C1, C2 := , ;

 > display(spacecurve(C1, x=-1..1), spacecurve(C2, x=-1..1),         color=[red,blue], thickness=5);

Surface that interpolates between the two space curves

 > C1 + <0, 0, z>; plot3d(%, z=C1[3]..C2[3], x=-1..1, lightmodel=light4);

PS: In this particular example the interpolating surface is vertical.  You may construct an interpolating surface even when the two curves are not situated one directly above the other. Then the interpolating surface will be:

```S := (1-z)*C1 + z*C2;
plot3d(S, z=0..1, x=-1..1, lightmodel=light4);```

## Here is how...

You didn't say maximum value of what, so here I show how to calculate the maxima (and minima) of S(t).  Adjust as needed.

Continuation of the code in  deplot3d-animated_SIA.mw.

Finding the critical points of S(t)

Here we will find the maxima (and minima) of S(t) when the solution starts at ICs8[1].

Change as needed.

The critical points of S occur when the derivative is zero.  So introduce a fourth

differential equation, dS/dt = U, and solve the system of four differential equations:

 > de4 := diff(S(t),t) = U(t);

 > IC := ICs8[1][];

 > dsol := dsolve({de1,de2,de3,de4, IC}, numeric, output=operator);

 > my_S := subs(dsol, S); my_U := subs(dsol, U);

 > plot([my_S(t), my_U(t)], t=0..6, color=["Red", "Green"]);

The times when S(t) has a local maximum or minimum:

 > t_vals := [fsolve(my_U(t), t=0..6, maxsols=10)];

 > S_vals := map(my_S, t_vals);

You may determine which is min and which is max visually.

It's possible to automate that by intorducing a fifth differential

equation dU/dt = V, and picking min and max based on the

sign of V at the cricial points, but perhaps that's more that what

you would care to do.

 >

deplot3d-animated_SIA-extra.mw

## The problem is ill-posed...

Your initial value problem is ill-posed.  Applying Dirac's delta at time

and specifying  at the same time is not quite meaningful.

To see why, apply Dirac's delta not at , but at some unspecified .

 > restart;
 > de := D(y)(t) + y(t)/tau = Dirac(t-a)/tau;

 > ic := y(0) = 0;

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

You are interested in the case of .  But the limit of the solution

as  goes to zero does not exist:

 > limit(sol, a=0, right);

 > limit(sol, a=0, left);

It seems that you favor one of the limits over the other but that

bias is only in your head;  the mathematics that you are applying

is not aware of it.

## Symbolic solution...

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.

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

 >

## A numerical solution...

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

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];  # . 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);