## 7776 Reputation

17 years, 361 days

## Without the goto...

I don't know how the undocumented goto works.  Here is a  version without goto.  Have a look at the help page for break.

 > restart;
 > for i to 4 do         for j to 4 do                 print([i,j]):                 if i = 2 and j = 3 then                         break i;  # break out of the 'i' loop                 end if;         end do; end do;

## Maybe this?...

I see a lot of experiments in your worksheet but I didn't see a clear statement of the objective.  Perhaps this is what you have in mind?

 > restart;
 > de := diff(phi(s), s, s) + K*cos(phi(s)) = 0;

 > bc := D(phi)(L) = phi__0;

 > dsol := dsolve({de,bc}, phi(s), implicit);

 > diff(dsol[1], s):  # either dsol[1] or dsol[2] will do isolate(%, diff(phi(s), s))^2;

## Here is how...

 > restart;
 > f := t -> sin(t) - 1/2;

 > plot(f(t), t=0..Pi);

 > plot([1, f(t) + 1], t=0..Pi, coords=polar, scaling=constrained, color=[blue,red], thickness=5);

## Do it in cylindrical...

What you have is not an error, but a bad choice of coordinates.  The intersection of a cone and plane is best described in cylindrical coordinates.  You will get a perfect ellipse that way.

You have the cone z=x^2+y^2, and the plane z = 3 + y/3.  Switch to cylindrical coordinates with x=r*cos(t), y=r*sin(t).   Then the equations of the cone and the plane become z = r^2 and z=3+1/3*r*sin(t).  At the intersection we have r^2 = 2 + 1/3*r*sin(t).  Solve this quadratic for r in terms of t.  Pick either of the two solutions. Call it r(t). Then the intersection curve is the ellipse described as a space curve:

[  r(t)*cos(t), r(t)*sin(t), 3 + 1/3*r(t)*sin(t)] ,    -Pi < t < Pi

## Perhaps this?...

Double-check your equation.  Perhaps "y" should be "u"?

If so, your boundary value problem does have a solution.  Here it is.

 > restart;
 > de := diff(u(x),x,x) + 1/x*diff(u(x),x) + u(x) = 6 - 9*x + x^2 - x^3;

Solve with the boundary condition .

We will deal with  afterward.

 > dsol := dsolve({de, u(1)=0});

We pick  by trial-and-error (there ought to

be a better way!) to enforce :

 > subs(c__2=-2.5948, dsol); plot(rhs(%), x=0..1, color=red, thickness=5);

## Maybe this?...

I can't see a meaning in what you have written.  Perhaps this is what you really want?

 > restart;
 > n := [ 82, 79, 84, 85, 76, 76, 78, 82, 85, 77, 83, 84, 87, 86, 87, 88, 76 ];

 > nops(n);

 > for i from 3 to nops(n) do         f[i] := n[i-2] - 3*n[i-1] + 3*n[i]; end do;

## Solution...

In your worksheet you have things like d*u/dt.  That says "d" times "u" divided by "dt". That's not how one expresses a derivative in Maple. Considering that, I suggest that you begin with learning about the Maple syntax, solve a few simple ordinary differential equations, and then a few simple partial differential equations, and only then attempt to solve a nonlinear system of PDEs.

For whatever it's worth, I have converted your sketches of equations to proper Maple syntax and solved them. I can't tell how useful will this be to you.

 > restart;
 > pde1 := diff(u(y,t),t) - s*diff(u(y,t),y) = diff(u(y,t),y,y) + delta*v(y,t);

 > pde2 := diff(v(y,t),t) - s*diff(v(y,t),y) = Pr*diff(v(y,t),y,y) + lambda*exp(v(y,t)/(epsilon*v(y,t)+1))/Pr;

 > ic := u(y,0)=0, v(y,0)=0;

 > bc := u(0,t)=0, v(0,t)=0, u(1,t)=1, v(1,t)=0;

 > s := 0.5; epsilon := 0.01; Pr := 7; lambda := 1; delta := 1;

 >
 > pdsol := pdsolve({pde1,pde2}, {ic, bc}, numeric, spacestep=0.005);

 > pdsol:-plot3d(u, t=0..0.3);

 > pdsol:-plot3d(v, t=0..0.3);

 >

## Numerical solution...

The second form of your differential equation has explicitly defined numerical coefficients except for V which is left undefined.  So we apply dsolve() to solve the differential equation numerically while leaving V as a parameters that can be assigned values retroactively.

 > restart;
 > de := diff(f(t), t) = (.21*.2894606222-.21*(1.158611162/((1+0.4836862955e-1*t)*(1/(1+0.4836862955e-1*t))^.1986136226))-.2894606222*f(t))/V;

 > ic := f(0) = .21;

 > dsol := dsolve({de,ic}, numeric, parameters=[V], output=operator);

 > F := eval(f, dsol);

 > dsol(parameters=[1]);   # set V = 1

 > F(4);

 > plot(F(t), t=0..20);

 > dsol(parameters=[2]);   # set V = 2

 > F(4);

 > plot(F(t), t=0..20);

 >

## Try this...

```for i from 1 to 24 do
print(i);
if irem(i,5) = 0 then
print("press enter to continue");
end if;
end do;```

## Correction...

That's because in the definition of F you have

`z=-2*cos(q)-2.5 .. -2.5-2*cos(q)`

That should be

```z=-2*cos(q)-2.5 ..  2.5-2*cos(q)
```

## Simplified version...

C_R has pointed out the source of the error.  You are attempting to calculate "le centre" of a paraboloid.  But a paraboloid has no "centre".

It seems to me that you are making things much more complicated than need be,  and plotting planes and quadratic surfaces through implicitplot3d is a poor choice.  These can be much better plotted as parametric surfaces through plot3d.  Here, for instance, is the paraboloid you wish to plot:

```L := 4:  K := 1.3*L:
display(
plot3d(< 2*s*cos(t), 3*s*sin(t), s^2 >,
s=0..2, t=-Pi..Pi, color="Orange"),
plot3d([s,t,0], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
plot3d([0,s,t], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
plot3d([t,0,s], s=-L..L, t=-L..L, transparency=0.6, color="DarkGray"),
arrow([0,0,0], [K,0,0], color="Red"),
arrow([0,0,0], [0,K,0], color="Green"),
arrow([0,0,0], [0,0,K], color="Blue")
);```

## Multispan beam...

This website is broken again and won't allow showing the contents of a worksheet.

## Streamlines...

Did you really mean that velocity field?  That velocity field is time-dependent.  Plotting time-dependent streamlines calls for an animation; a static plot won't do.   Usually when people talk about streamlines, they mean the streamlines of a stationary (that is, time-independent) velocity field.

Here I plot the streamlines of a stationary velocity field.  This may be adapted to produce an animation of the time-dependent case if necessary, but that would call for a closer familiarity with Maple.

 > restart;
 > with(DEtools):

Velocity field:

 > v := (x,y) -> <1+x^2+y^2, x+y>;

Streamlines are the solutions of this system of differential equations

 > de := [seq](diff(,t) =~ v(x(t),y(t)));

Here is the corresponding vector field:

 > DEplot(de, [x(t),y(t)], t=0..1, x=-3..3, y=-3..3);

Specify initial conditions along the left edge:

 > ic := seq([x(0)=-3, y(0)=h], h=-3..3, 0.2);

Plot streamlines

 > DEplot(de, [x(t),y(t)], t=0..3, [ic], x=-3..3, y=-3..3, linecolor=black, thickness=1);

## Rotating helix...

This worksheet provides an alternative approach to what I had provided earlier in RevolSurfaceQ.mw and adds animation as you have requested.

 > restart;
 > with(plots):

rotates the vector  about the vector  through an angle ,

The length of the vector  is immaterial as long as it is nonzero.

 > R := proc(u::Vector(3))         local pname, n, phi, nlen, t;         pname := op(procname);         if not type(procname, 'indexed') or nops([pname])<>2 then                 error "Two indices needed";         end if;         n := pname[1];         phi := pname[2];         if not type(n, 'Vector(3)') then                 error "expected Vector(3) for the first argument of R[n,phi]";         end if;         nlen := simplify(sqrt(n^+ . n));         if nlen = 0 then                 error "the vector n in R[n,phi] should have a nonzero length";         end if;         n /= nlen;         t := (n^+ . u)*n;         t + (u - t)*cos(phi) + LinearAlgebra:-CrossProduct(n,u)*sin(phi); end proc:

Example

Let's rotate the helix

 > H := < cos(s), sin(s), s/4 >;

 > L := t -> ;

 > n := <1, 0, 1>;

through an angle .  We get

 > S := R[<1,0,1>,q](H);

For a fixed value of  and varying , the vector  calculated above
represents the parametric equation of the rotated helix.

If both  and  vary, the vector  represents the parametric equation

of the surface swept by the helix as it rotates about the vector .

Here is what the rotating helix looks like:

 > demo1 := q -> display(         tubeplot([seq](R[n,q](H)), s=-4*Pi..4*Pi, radius=0.08, numpoints=200, color="Green"),         tubeplot([seq](L(s)), s=-2..2.5, radius=0.08, color="Red")):
 > animate(demo1, [q], q=0..2*Pi, lightmodel=light4, style=surface,         orientation=[-120,75,0], labels=[x,y,z]);

and here is the rotating helix and the surface it sweeps:

 > S := a -> display(         plot3d(R[n,q](H), s=-4*Pi..4*Pi, q=0..a, color="Orange"),         tubeplot([seq](R[n,a](H)), s=-4*Pi..4*Pi, radius=0.08, numpoints=200, color="Green"),         tubeplot([seq](L(s)), s=-2..2.5, radius=0.08, color="Red") ):
 > animate(S, [a], a=0..2*Pi, style=surface, lightmodel=light4,         orientation=[-120,75,0], labels=[x,y,z], frames=60, axes=framed);

## Animation of the ruled surface...

Here is an animation of the generatrix of your ruled surface.

And here is the worksheet the produced it.