## 13776 Reputation

14 years, 64 days

## I'm going to take a wild guess...

that the attached is what the OP was aiming for

 >
 > restart;   _EnvHorizomtalName='x':   _EnvVerticalName='y':   aProc:=proc(t)               uses geometry, plots:               local a := 11,                     b := 7,                     R := sqrt(a^2 + b^2),                     lt:= y=m*x+c,                     sol, l1,l2, ell, xv, yv;               point( P, [R*sin(t), R*cos(t)]);               point( OO, [0,0]);               circle( cir, [OO, R]);               ellipse(ell, x^2/a^2+y^2/b^2-1=0, [x, y]):               sol:= [ solve                       ( subs                         ( c=VerticalCoord(P)-m*HorizontalCoord(P),                           discrim                           ( lhs                             ( collect                               ( expand                                 ( subs                                   ( lt,                                     Equation(ell)                                   )                                 ),                                 x                               )                             ),                             x                           )                         )=0,                         m                       )                     ];               line( l1,                     eval                     ( y=m*x-m*HorizontalCoord(P)+VerticalCoord(P),                       m=sol[1]                     ),                     [x,y]                   );               line( l2,                     eval                     ( y=m*x-m*HorizontalCoord(P)+VerticalCoord(P),                       m=sol[2]                     ),                     [x,y]                   );               xv:=evalf~                   ( evalc~                     ( [ solve                         ( subs                           ( isolate(Equation(l1),y),                             Equation(ell)                           )                         )                       ]                     )                   )[1];               yv:= solve                    ( eval                      ( Equation(l1),                        x=xv                      )                    ) ;               point( S, [xv, yv]);               xv:=evalf~                   ( evalc~                     ( [ solve                         ( subs                           ( isolate(Equation(l2),y),                             Equation(ell)                           )                         )                       ]                     )                   )[1];               yv:= solve                    ( eval                      ( Equation(l2),                        x=xv                      )                    );               point( T, [xv, yv]);               display               ( [ textplot                   ( [ [ coordinates(P)[], "P"],                       [ coordinates(S)[], "S"],                       [ coordinates(T)[], "T"]                     ],                     font=[times, bold, 16],                     align=[above, right]                   ),                                           draw                   ( [ ell(color=blue),                       cir(color=blue),                       l1(color=green),                       l2(color=red),                        S(color=black, symbol=solidcircle, symbolsize=16),                        T(color=black, symbol=solidcircle, symbolsize=16),                        P(color=black, symbol=solidcircle, symbolsize=16)                     ],                     scaling=constrained,                     axes=none,                     view=[-15..15, -15..15]                   )                 ]               );         end proc:
 > nFig := 60:  plots:-display(seq(aProc(2*Pi*i/nFig), i=0..nFig), insequence = true);
 >

## I still...

don't really understand what you are trying to achieve, or why, but maybe something in the code below will help

```  restart;
with(GraphTheory):

L1:=[ [ (0, 1), (1, 2), (1, 10), (2, 3), (3, 4), (4, 5), (4, 9), (5, 6), (6, 7), (7, 8),
(8, 9), (10, 11), (11, 12), (11, 16), (12, 13), (13, 14), (14, 15), (15, 16)
],
[ (0, 10), (1, 2), (1, 9), (2, 3), (3, 4), (4, 5), (4, 9), (5, 6), (6, 7), (7, 8),
(8, 9), (10, 11), (11, 12), (11, 16), (12, 13), (13, 14), (14, 15), (15, 16)
]
]:
g:= [ seq
( Graph
( { seq
( {L1[j][i-1], L1[j][i]},
i=2..numelems(L1[j]),2
)
}
),
j=1..numelems(L1)
)
]:

g[1];
Edges(g[1]);
g[2];
Edges(g[2]);
DrawGraph(g[1]);
DrawGraph(g[2]);
```

I would normally upload a worksheet , but someone seems to have broken the big green uparrow. Let's hope it is (very) temporary

## Something(?)...

The first thing to realise is that pdsolve(...., numeric) does not return any information about the derivatives of the dependent function T(y,t). The only way I can think of to generate this information is to use the pds:-value() method combined with the numeric differentiation command fdiff(), to obtain derivative values at a point.

Unfortunately the fdiff() command is very slow, so in the attached, I have restricted the computation of your desired functions

1+Nr*(T(y, t)+1)^3)*(diff(T(y, t), y)), at t = 1,

diff(T(y, t), y) at y = 0

diff(T(y, t), y) at y = 1

to only a few points across the range of the relevant variable. Even with this restriction, the attached worksheet takes several minutes to execute on my machine.

The last plot in the attached, ie diff(T(y, t), y) at y = 1 for Nr=[0, 1, 10] essentially produces zero for all t and any Nr. I am not sure whether this is plausible or not!

 > restart;    with(plots):   PDE1 := Pr*(diff(T(y, t), t)-Ree*(diff(T(y, t), y))) = (1+Nr*(T(y, t)+1)^3)*(diff(T(y, t), y, y))+3*Nr*(T(y, t)+1)^2*(diff(T(y, t), y))^2;   ICandBC := {T(1, t) = 1, T(y, 0) = 1, (D[1](T))(0, t) = T(0, t)}:   Ree := .1:   Pr := 6.2:   HA1 := [0, 1, 10]:   AA := [red, green, blue, cyan, purple, black]:   for i to nops(HA1) do      Nr := op(i, HA1);      print("Nr = ", %);      PDE[i] := {PDE1};      pds[i] := pdsolve( PDE[i],                         ICandBC,                         numeric,                         spacestep = 1/200,                         timestep = 1/100                       );      PlotsT[i] := pds[i]:-plot( T(y, t),                                 t = 1,                                 linestyle = "solid",                                 labels = ["y", "u"],                                 color = AA[i],                                 numpoints = 800                               );      Tv:=eval( T(y,t),                pds[i]:-value( t=1,                               output=listprocedure                             )              ):      f:=yv->(1+Nr*(Tv(yv)+1)^3)*fdiff(Tv(y), y=yv):      plt1[i]:=plot( 'f(y)',                    y=0..1,                    color=AA[i],                    size = [1000, 600],                    axes = boxed,                    axesfont = ["Arial", 14, Bold],                    thickness = 3                  ):      Tv:=tv-> eval(T(y,t), pds[i]:-value( y=0, t=0..1 )(tv));      f:=tau->fdiff('Tv(t)', t=tau);      plt2[i]:= plot( 'f(t)',                      sample=[seq(j, j=0..1, 0.1)],                      adaptive=false,                             color=AA[i],                      style=point,                      view=[0..1, default],                      symbol=solidcircle,                      symbolsize=16,                      size = [1000, 600],                      axes = boxed,                      axesfont = ["Arial", 14, Bold],                      thickness = 3                   );      Tv:=tv-> eval(T(y,t), pds[i]:-value( y=1, t=0..1 )(tv));      f:=tau->fdiff('Tv(t)', t=tau);      plt3[i]:= plot( 'f(t)',                      sample=[seq(j, j=0..1, 0.1)],                      adaptive=false,                             color=AA[i],                      style=point,                      view=[0..1, default],                      symbol=solidcircle,                      symbolsize=16,                      size = [1000, 600],                      axes = boxed,                      axesfont = ["Arial", 14, Bold],                      thickness = 3                   );   od:      display( [seq(PlotsT[k], k = 1 .. nops(HA1))],            size = [1000, 600],            axes = boxed,            labels = [x, (convert("T", symbol))(x, T)],            labelfont = ["Times", 14, Bold],            labeldirections = [horizontal, vertical],            axesfont = ["Arial", 14, Bold],            thickness = 3          );   display( [seq(plt1[j], j=1..nops(HA1))],            title=typeset((1+'Nr'*(T(y, t)+1)^3)*(diff(T(y, t), y)), "at t = 1"),            titlefont=["Arial", 16, Bold]          );   display( [seq(plt2[j], j=1..nops(HA1))],            title=typeset( diff(T(y, t), y), "at y = 0"),            titlefont=["Arial", 16, Bold]          );   display( [seq(plt3[j], j=1..nops(HA1))],            title=typeset( diff(T(y, t), y), "at y = 1"),            titlefont=["Arial", 16, Bold]          );
 >

## Leaving aside...

minor syntax issues, the big problem you have is that your PDE is

1. second-order in 'x' ,requiring two boundary conditions, and
2. first-order in 't' requiring one initial condition.

You only have a total of two conditions

## Not a "Hard-coded number"...

when you write X[`1`], the backquotes convert whatever is between them to a name., so you now have a variable whose name is 1, and whose value can be anything you want. You can tnen set the variable named `1` to be the value 1, as shown in the attached.

By the way this only "works" if you use the long-deprecated (since Maple 8 in 20002) vector() constructor. If you had used the recommended Vector() constructor, then X[`1`] would generate an error.

 > restart;   X := vector([seq(x, x = 1 .. 5)]);   Y := vector([seq(x, x = 6 .. 10)]);   X[1];   Y[1];   Yfactor := evalf(X[`1`])/Y[`1`];   eval(Yfactor, `1`=1);   eval(Yfactor, `1`=2);   eval(Yfactor, `1`=3);
 (1)
 >

## Use backquotes...

See the "toy" example in the attached

 > plot( x^2,       x=0..1,       labels=[ typeset("a-value=", `A__*`),                typeset("y-value=", `A__*`^2)              ]     );
 >

## I suspect...

that your basic problem is that eqn3 contains the parameter 'T' which is defined nowhere. In the attached I have just set this equal to 1.

There are various other "inconsistencies" in the code, eg use of non-commutative multiplication, use of greek symbols together with their "roman" transliteration. Theses latter two are probably(?) not significant, but since I wouln't trust MAple's 2D input parser further than I could throw it, I fixed these up as well

At least the attached now executes!

 >
 >
 >
 >
 >
 >
 (1)
 >
 (2)
 >
 (3)
 >
 (4)
 >
 (5)
 >
 (6)
 >
 (7)
 >
 (8)
 >
 (9)
 >
 (10)
 >
 (11)
 >
 (12)
 >
 >
 (13)
 >
 >
 >
 >
 >
 >

## Have you considered...

the Maple command ssystem(command) which passes command to the host operating system which, in turn, performs the appropriate function?

See the help at ?ssystem

## Physics:-Diff()...

is *probably* the best way, but you could aslo jut "cheat" a little. Substitute theta(t)=theta in you expression, differentiate with respect to theta, and the (if necessary/desired) substitute theta=thta(t) in the result

As in the attached

 > restart;   p := (m1*a1 + m2*(a123 + z(t)))*sin(theta(t))*g + m3*(sin(theta(t))*(a1234 + z(t)) + cos(theta(t))*a5)*g: # # Either #   diff   ( eval     ( p,       theta(t)=theta     ),     theta   ); # # or, #   eval   ( diff     ( eval       ( p,         theta(t)=theta       ),       theta     ),     theta=theta(t)   );
 (1)
 >

## Well, some simple typos...

for example, in the draw() command , you want to draw

seg1,seg2,seg3,seq4,

there is a difference between seg4 and seq4". Also there is no command

Triangle()

although there is a command

triangle()

Still on the subject of simple typos, you might want to consider the spelling of

_EnvHorizont:lName = 'x';

I fixed all of these in the attached, and probably a few more things. I also deleted several lines of completely redundant coded

 > restart:   with(geometry):   with(plots):   _EnvHorizontalName = 'x':   _EnvVerticalName = 'y':    R := 5:    ang := [3/4*Pi, -(3*Pi)/4, -Pi/6,4*Pi/9]:    seq    ( point      ( `||`(P, i),        [ R*cos(ang[i]), R*sin(ang[i])]      ),      i = 1 .. 4    ):    seq    ( dsegment      ( `||`(seg, i),        [ `||`(P, i),          `||`(P, irem(i, 4) + 1)        ]      ),      i = 1 .. 4    ):    triangle(Tr1,[P1,P2,P4]):    EulerCircle(Elc1,Tr1,'centername'=o):    circle(cir, [point(OO, [0, 0]), R]):    display    ( [ draw        ( [ P1(color = black, symbol = solidcircle, symbolsize = 12),            P2(color = black, symbol = solidcircle, symbolsize = 12),            P3(color = black, symbol = solidcircle, symbolsize = 12),            P4(color = black, symbol = solidcircle, symbolsize = 12),            seg1,            seg2,            seg3,            seg4,            Tr1(color=green),            Elc1,            cir(color = blue)          ]        ),        textplot        ( [ seq            ( [ coordinates(`||`(P, i))[],                convert(`||`(P, i), string)              ],              i=1..4            )          ],          align=[above, right]        )      ],      axes=none    );
 >

## Mugged by 2D input...

Although you do not post a worksheet (again), I can tell that you are using Maple's 2-D input mode. You have to be careful when doing this, because whitespace is *sometimes* interpreted as an implicit multiplication. However you managed to post code here (some kind of cut-and-paste, I guess) actually illustrates the problem

```display*([draw*[P1(color = black, symbol = solidcircle, symbolsize = 12),
P2(color = black, symbol = solidcircle, symbolsize = 12),
P3(color = black, symbol = solidcircle, symbolsize = 12),
cir(color = blue)],
textplot*([[coordinates(P1)[], "P1"],
[coordinates(P2)[], "P2"],
[coordinates(P3)[], "P3"]], align = [above, right])], axes = none);```

Notice that instead of display(, you have display*( And the same with the 'draw and 'textplot' commands.

This problem is fixed in the attached which uses Maple's 1D-input mode where whitespace is whitespace and thus the user can use as much indentation (ie spaces) as desired, in order to make code easily readable. I also deleted some code whihc you don't seem to be using.

 > restart;   with(geometry):   with(plots):   _EnvHorizontalName = 'x':   _EnvVerticalName = 'y':   R := 5:   ang := [2/3*Pi, -3*Pi*1/4, -Pi*1/6]:   seq   ( point     ( `||`(P, i),       [ R*cos(ang[i]), R*sin(ang[i]) ]     ),     i = 1 .. 3   ):   seq   ( dsegment     ( `||`(seg, i),       [ `||`(P, i), `||`(P, irem(i, 3) + 1) ]     ),     i = 1 .. 3   ):   circle   ( cir,     [ point(OO, [0, 0]), R ]   ):   display   ( [ draw       ( [ P1(color = black, symbol = solidcircle, symbolsize = 12),           P2(color = black, symbol = solidcircle, symbolsize = 12),           P3(color = black, symbol = solidcircle, symbolsize = 12),           cir(color = blue)         ]       ),       textplot       ( [ [coordinates(P1)[], "P1"],           [coordinates(P2)[], "P2"],           [coordinates(P3)[], "P3"]         ],         align = [above, right]       )     ],     axes = none   );
 >

## Your basic problem...

is that you have

``` x := C[1] + n[1]*t;
y := C[2] + n[2]*t;
EQ := eq = 0;
tt := solve(subs(x = C[1] + n[1]*t, y = C[2] + n[2]*t, EQ), t);```

Since 'x' and 'y' have been assigned values,  the subs() commmand will actually evaluate to

`subs( C[1] + n[1]*t= C[1] + n[1]*t,  C[2] + n[2]*t= C[2] + n[2]*t, EQ)`

which I'm guessing isn't what you want

Maybe you intended something like the attached?

PS if you really want to know the projection of a point on to a line, you can do this with the geometry:-projection() command

 > restart; with(LinearAlgebra): A := [1, -2]: B := [-2, 3]: C := [1, 1]: M := [x, y]: ProjPL := proc(C, A, B)        local M, AB, AM, Q, eq, EQ, eq1, a, b, c, t, tt, n, dist,              x, y, xH, yH, H, CH, no;              M := [x, y];              AM := M - A;              AB := B - A;              Q := Matrix(2, [AM, AB]);              eq := Determinant(Q);              a := coeff(eq, x);              b := coeff(eq, y);              c := tcoeff(eq);              dist := abs(a + b + c)/sqrt(a^2 + b^2);              n := [a, b];              EQ:=eq = 0;              tt := solve(subs([x=C[1] + n[1]*t,y=C[2] + n[2]*t], EQ), t);              x := C[1] + n[1]*t;              y := C[2] + n[2]*t;              xH := subs(t = tt, x);              yH := subs(t = tt, y);              H := [xH, yH];              CH := H - C;              no := sqrt(CH[1]^2 + CH[2]^2);              return EQ, dist, H, no;         end proc: ProjPL(C, A, B);
 (1)
 >

## Anything wrong...

with the most basic use of the dsolve() command, as in the attached

 > restart:   odeSys:= diff(y__1(x),x)=y__3(x)-cos(x),            diff(y__2(x),x)=y__3(x)-exp(x),            diff(y__3(x), x)=y__1(x)-y__2(x):   ics:= y__1(0)=0,         y__2(0)=0,         y__3(0)=0:   dsolve([odeSys, ics]);
 (1)
 >

## To get the desired form...

is a bit more complicated - see the attached

 > restart;   u(x,t):=;   pde:=diff(u(x,t), t) = Matrix([[0, mu*k^2/2], [2*A^2 - mu*k^2/2, 0]]).u(x,t);   ans:=pdsolve        ( simplify~          ( subs            ( isolate              ( w^2 = -mu^2*k^4/4 + mu*k^2*A^2,                mu               ),               pde            )          )        );
 (1)
 >

## Still no worksheet!!!...

In the absence of a worksheet, I have to hack your cut+paste entry to get something executable. It generates the error message

Error, (in geometry:-triangle) The three given points are AreCollinear

This error message is produced by the command

triangle(TR, [M1, M2, A]):

because two of these points (A and M1) have the same coordinates, so no triangle can be constructed. Coordinates of all three points are shown below

A=[2, 4.582575695]
M1=[2.000000000, 4.582575695]
M2=[1.330063498, -4.819847621]

Since your algorithm is therefore incorrect, and I don't know what the "poncelet theorem for the triangle" is, the error is impossible for me to fix It *seems* as if you are doing something vaguely related to Poncelet's Closure Theorem (aka Poncelet's porism) described in

https://en.wikipedia.org/wiki/Poncelet%27s_closure_theorem

so, just for amusement I have included an illustration of this Theorem in the worksheet attached. It *may* be of use in whatever construction you are trying to produce

 > restart;   with(geometry):   with(plots):   _EnvHorizontalName='x':   _EnvVerticalName='y':   R := 5.0:   ang := evalf~([2/3*Pi, -3*Pi*1/4, -Pi*1/6]):   seq( point( `||`(P, i),               [ R*cos(ang[i]),                 R*sin(ang[i])               ]             ),        i = 1 .. 3      ):   triangle( T1,             [seq(`||`(P, j), j=1..3)]           ):   incircle(C1, T1):   circumcircle(C2, T1):   porism:= proc(t)                 local theta:=ang[1]+t,                       P1, P2, P3, L, II, T1,j;               #               # local point P1 is somewhere on the (global)               # circumcircle C2               #                 point( P1,                        [ R*cos(theta),                          R*sin(theta)                        ]                      ):               #               # Get the tangent lines from (local) P1 to               # the (global) incircle C1               #                 TangentLine(L, P1, C1):               #               # Find the intersection of these tangent lines               # with the (global) circumcircle C2. Note that               # each tangent line, intersects C2 at two points,               # one of which is (local) point(P1). So for               # both tangent lines, determine the intersection               # point whihc isn't P1               #                 intersection(II, L[1], C2):                 if   distance(II[1], P1) > distance(II[2], P1)                 then point( P2,                             coordinates(II[1])                           );                 else point( P2,                             coordinates(II[2])                           );                 fi:                 intersection( 'II', L[2], C2):                 if   distance(II[1], P1) > distance(II[2], P1)                 then point( P3,                             coordinates(II[1])                           );                 else point( P3,                             coordinates(II[2])                           );                 fi:                 triangle( T1,                           [P1, P2, P3]                         ):                 display                 ( [ draw                     ( [ P1(color=black, symbol=solidcircle, symbolsize = 12),                         P2(color=black, symbol=solidcircle, symbolsize = 12),                         P3(color=black, symbol=solidcircle, symbolsize = 12),                         T1(color=red),                         C1(color=blue),                         C2(color=green)                       ]                     ),                     textplot                     ( [ [ coordinates(P1)[], "P1"],                         [ coordinates(P2)[], "P2"],                         [ coordinates(P3)[], "P3"]                       ],                       align=[above, right]                     )                   ],                   axes=none                 );         end proc:   nf:=60:   display   ( [ seq       ( porism(j),         j=0..evalf(2*Pi*(1-1/nf)), evalf(2*Pi/nf)       )     ],     insequence=true   );
 >