## 4993 Reputation

9 years, 213 days

## With minimal changes...

to your original code, one can produce the attached

 > with(plottools): with(plots): display( line([.8, 0], [1, .2], color = red),          line([.6, 0], [1, .4], color = red),          line([.4, 0], [1, .6], color = red),          line([.2, 0], [1, .8], color = red),          line([ 0, 0], [1,  1], color = red),          line([0, .2], [.8, 1], color = red),          line([0, .4], [.6, 1], color = red),          line([0, .6], [.4, 1], color = red),          line([0, .8], [.2, 1], color = red),          rectangle([0, 1], [1, 0], color=white, thickness = 1),          axes=none        ) >

## Is this...

what you are trying to achieve???

By the way this should be a "Question" not a "Post"

 > restart;   with(plots):   sol:=solve( [ x^2-y=0, y-(x+5)=0],               [x, y],               explicit             ):   display( [ inequal( [ x^2-y<0,                         y-(x+5)<0                       ],                       x=-4..4,                       y=0..10,                       color=red                     ),              textplot( [ [ rhs~(sol)[],                            "A",                            align=right                          ],                          [ rhs~(sol)[],                            "B",                            align=left                          ]                        ],                        font=[times, bold, 18]                      )            ],            size=[600, 600]          ); > ## Use 'singsol' option...

as in the attached. See the help at ?dsolve/details

I don't know why it is sometimes necessary to use the 'singsol' option, and sometimes the singular solution appears without being specifically requested. However if you are partcularly interested in singular solutions, then it is probably a good idea to make this requirement explicit.

 > restart;
 > Typesetting:-Settings(typesetprime=true):
 > ode:=x^2*diff(y(x),x)^2-(1+2*x*y(x))*diff(y(x),x)+1+y(x)^2 = 0; (1) (2)
 > Vector([dsolve(ode,y(x))]); #now it shows singular solution (first one below) (3)
 > PDEtools:-casesplit(ode) (4)
 > ode:=convert(ode,D): #solve for y(x) first, this will generate 2 ODE's sol:=[solve(ode,y(x))]: odes:=Vector(map(z->y(x)=z,convert(sol,diff))) (5)  (6)
 > dsolve(odes,y(x)); #where is singular solution? dsolve(odes,y(x), singsol=all);  #here  (7)
 > dsolve(odes,y(x)); #where is singular solution? dsolve(odes,y(x), singsol=all);  #here  (8)
 >

## The easy way...

is just to use Maple's builtin numeric solvers.

The attached reproduces the graphs in Figure 1 and Figure 2 of the reference you supply.

 > restart;   with(plots): # # Set up odes and BCs #   odes:=(1+K)*diff(u(y),y,y)+K*diff(N(y),y)+theta(y)=M*u(y),         (1+K/2)*diff(N(y),y,y)-K*(2*N(y)+diff(u(y),y))=0,         diff(theta(y),y,y)=0:   bcs:=u(0)=0, theta(0)=R, N(0)=0,        u(1)=0, theta(1)=R, N(1)=0:
 > # # Solve system for various values of the parameter 'M' #   mVals:=[0,5,10]:   cols:=[red, blue, green]:   solM:=[ seq           ( dsolve             ( eval               ( [odes, bcs],                 [K=5, R=0.5, M=i]               ),               numeric             ),             i in mVals           )         ]:
 > # # Display u(y) versus y and N(y) versus y for # each of the values of M #   display( [ seq              ( odeplot                ( solM[i],                  [y, u(y)],                  y=0..1,                  color=cols[i],                  legend=typeset("M=",mVals[i]),                  legendstyle=[font=[times, bold, 16]]                ),                i=1..numelems(solM)              )            ],            axes=boxed,            title=typeset( "Effect of M on velocity profile (u)"),            titlefont=[times, bold, 20],            size=[800,800]          );   display( [ seq              ( odeplot                ( solM[i],                  [y, N(y)],                  y=0..1,                  color=cols[i],                  legend=typeset("M=",mVals[i]),                  legendstyle=[font=[times, bold, 16]]                ),                i=1..numelems(solM)              )            ],            axes=boxed,            title=typeset( "Effect of M on microrotation (N)"),            titlefont=[times, bold, 20],            size=[800,800]          );  > # # Solve system for various values of the parameter 'K' #   KVals:=[0,5,10]:   solK:=[ seq           ( dsolve             ( eval               ( [odes, bcs],                 [K=i, R=0.5,M=5]               ),               numeric             ),             i in KVals           )         ]:
 > # # Display u(y) versus y and N(y) versus y for # each of the values of K #   display( [ seq              ( odeplot                ( solK[i],                  [y, u(y)],                  y=0..1,                  color=cols[i],                  legend=typeset("K=",KVals[i]),                  legendstyle=[font=[times, bold, 16]]                ),                i=1..numelems(solK)              )            ],            axes=boxed,            title=typeset( "Effect of K on velocity profile (u)"),            titlefont=[times, bold, 20],            size=[800,800]          );   display( [ seq              ( odeplot                ( solK[i],                  [y, N(y)],                  y=0..1,                  color=cols[i],                  legend=typeset("K=",KVals[i]),                  legendstyle=[font=[times, bold, 16]]                ),                i=1..numelems(solK)              )            ],            axes=boxed,            title=typeset( "Effect of K on microrotation (N)"),            titlefont=[times, bold, 20],            size=[800,800]          );  >

## Just realised...

you want the triangle ABC in the x-y plane, with the origin at the centre of its circumcircle. The attached worksheet will do this, and then compute the relevant tetrahedron.

The only "tricky" thing about this is that part of the calculation is done in in the (2D) geometry() package, and the remainder is done with the 3D geom3d() package

 > restart;
 > with(geometry):
 > # # Define the sidelengths #   sa, sb, sc, ab, bc, ac:= 3,5,7,3,4,5:
 > ####################################################### # Define three points which will become the base triangle. # # Without loss of generality one can specify # # 1. point A2D is at the origin # 2. point B2D is along the +ve x-axis # 3. point C2D is in the x-y plane, with positive #    y-coordinate   point( A2D,          [0, 0]        ):   point( B2D,          [x__b, 0]        ):   point( C2D,          [x__c, y__c]        ): # # Using the given side lengths, solve the distance # equations, to obtain (and assign) the unknown # coordinates x__b, x__c, y__c #   sol:=  assign          ( solve           ( [ x__b > 0,               y__c > 0,               distance(A2D,B2D)=ab,               distance(A2D,C2D)=ac,               distance(B2D,C2D)=bc             ]           )         ):   triangle(T1, [A2D,B2D,C2D]): # # Get the circumcircle for the triangle and generate # the directed line segment between the centre of this # circumcircle and the origin #   circumcircle(cc2d, T1):   dsegment( ds,             point             ( CS,               coordinates               ( center                 ( cc2d                 )               )             ),             point             ( O2,               [0,0]             )           ): # # Translate teh circumcircle so that the centre # is at the origin. Translate the triange by the # same amount #   translation( cc2d_trans,                cc2d,                ds              ):   triangle( T1_trans,             [ translation               ( A2D_trans,                 A2D,                 ds               ),               translation               ( B2D_trans,                 B2D,                 ds               ),               translation               ( C2D_trans,                 C2D,                 ds               )              ]           ):   draw( [T1_trans, cc2d_trans]); # # Convert the coordinates of the triangle in 2D # to 3D, just by adding a '0' for the z-coorcinate. # These coordinates will be used later # # Unload the geometry package, ready for loading # the geom3d package #   cA:=[coordinates(A2D_trans)[],0];   cB:=[coordinates(B2D_trans)[],0];   cC:=[coordinates(C2D_trans)[],0];   unwith(geometry):    (1)
 > with(geom3d): # # Define the points for the tetrahedron. Note that # the coordinates for the points A, B, C are obtained # by solving the 2D problem above #   point(A, cA):   point(B, cB):   point(C, cC):   point(S, [xs, ys, zs]): # # Using the supplied distances, solve for the coordinates # of the "unknown" point S #   sol:=  assign          ( solve           ( [ zs>0,               distance(S,A)=sa,               distance(S,B)=sb,               distance(S,C)=sc             ]           )         ): # # Define a general tetrahedron based on the points A, B, C, S # where by construction A,B,C will lie in the x-y plane, and # the origin will be the centre of the circumcircle of the # triangle ABC #    gtetrahedron(tet,                 [A,B,C,S]                ):    draw(tet, axes=boxed, scaling=constrained);    detail(tet);  (2)
 >

## Volume of a terahedron...

using the Cayley-Menger determinant and the geom3d() package

 > restart;
 > with(geom3d):   with(LinearAlgebra):
 > # # Define the sidelengths #   sa, sb, sc, ab, bc, ac:= 3,5,7,3,4,5: # # Compute the volume from the OP-supplied formula #   V1:= (1/12)*sqrt( sa^2*bc^2*(ab^2+ac^2-bc^2-sa^2+sb^2+sc^2)                     +                     sb^2*ac^2*(ab^2-ac^2+bc^2+sa^2-sb^2+sc^2)                     +                     sc^2*ab^2*(-ab^2+ac^2+bc^2+sa^2+sb^2-sc^2)                     -                     sa^2*ab^2*sb^2                     -                     sb^2*bc^2*sc^2                     -                     sa^2*ac^2*sc^2                     -                     ab^2*bc^2*ac^2                  ); (1)
 > ################################################################# # Construct the Cayley-Menger matrix for 4 points in 3-space # from the Wikipedia entry at # # https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant #   CM:= Matrix( 3,                3,                [ [ 2*d[0,1]^2,                 d[0,1]^2+d[0,2]^2-d[1,2]^2, d[0,1]^2+d[0,3]^2-d[1,3]^2 ],                  [ d[0,1]^2+d[0,2]^2-d[1,2]^2,           2*d[0,2]^2,       d[0,2]^2+d[0,3]^2-d[2,3]^2 ],                  [ d[0,1]^2+d[0,3]^2-d[1,3]^2, d[0,2]^2+d[0,3]^2-d[2,3]^2,          2*d[0,3]^2]                ]              ); # # Substitute the OP's side lengths in the above matrix and # compute the determinant, hence the enclosed volume. Note # that this volume must be positive for the four points # to define a "realisable" tetrahedron # # Note that the computed volume agrees with that found by # the OPs formula above and that returned by the geom3d() # package, later #   V2:= sqrt        ( Determinant          ( eval            ( CM,              [ d[0,1]=sa,                d[0,2]=sb,                d[0,3]=sc,                d[1,2]=ab,                d[1,3]=ac,                d[2,3]=bc              ]            )          )/288        );  (2)
 > ####################################################### # Define four points which will become the tetrahedron. # # Without loss of generality one can specify # # 1. point A is at the origin # 2. point B is along the +ve x-axis # 3. point C is in the x-y plane, with positive #    y-coordinate # 4. point S is in the positive z-hemisphere #   point( A,          [0, 0, 0]        ):   point( B,          [x__b, 0, 0]        ):   point( C,          [x__c, y__c, 0]        ):   point( S,          [x__s, y__s, z__s]        ): # # Using the given side lengths, solve the distance # equations, to obtain (and assign) the unknown # coordinates x__b, x__c, y__c, x__s, y__s, z__s #   sol:=  assign          ( solve           ( [ x__b > 0,               y__c > 0,               z__s > 0,               distance(A,B)=ab,               distance(A,C)=ac,               distance(B,C)=bc,               distance(S,A)=sa,               distance(S,B)=sb,               distance(S,C)=sc             ]           )         ): # # Create the the tetrahedron from the four points # and the sphere which goes through its vertices #   gtetrahedron( tet,                 [A, B, C, S]               ):   sphere( sph,           [A, B, C, S]         ): # # Return the volume of the tetrahedron. Note that this # agrees with that obtained from both the OP-supplied # formula and the Cayley-Menger determinant given above # # NB slightly surprised that this "worked", because # Maple documentation suggests that the volume() method # only works for "regular" polyhedra (ie objects of # type 'tetrahedron3d', but not objects of type # 'gtetrahedron3d') #   V3:=volume(tet); # # Get the vertices of the tetrahedron #   vertices(tet);  (3)
 > ################################################## # Generate the directed line segment between the # center of the sphere and the origin #   dsegment( ds,             point             ( CS,               coordinates               ( center                 ( sph                 )               )             ),             point             ( OO,               [0,0,0]             )           ): # # Translate the sphere found earlier so that its # center is at the origin. # # Perform the same translation on the tetrahedron #   translation( sph_trans,                sph,                ds              ):   translation( tet_trans,                tet,                ds              ): # # Draw the tranlated sphere and tetrahedron #   draw( [ sph_trans           ( style=surface,             color=red,             transparency=0.9           ),           tet_trans           ( style=surface,             color=black           )         ],         axes=boxed       ); # # Get the vertices of the translated tetrahedron #   vertices(tet_trans);  (4)
 >

## Not sure...

that a 3D-plot would be a good idea - anything wrong with producing a 2-D plot for each differetn value of the variable 'tswitch' as in the attached? Set up functions

 >  Start plotting

 >   >  (2.1)
 > ## Some extensions...

to Rouben's solution, which incorporate

1. Using the Optimization package to obtain the desired diffusion coefficient: produces the same value as Rouben, ie 0.138865421347717
2. Separating the setup and solving of the PDE from the routine used to fit it. This makes it easier to reuse the same code to generate the complete solution when the desired diffusion coefficient is found from (1) above
3. From (2) above, one can calculate the residual errors for each of the values in the vector t_number, as well as producing plots etc, using the solution with the optimal diffusion coefficient found in (1) above

See the attached (based on Rouben's code)

 > restart;
 > # # Parameters and procedures #   t_number:=<0, 10, 20, 30, 40>:   m_number:=<11.50000000, 4.641182547, 1.273311993, 0.361845238, 0.288711649>:   L:=2:   getSol:= proc(d)                 local Mx0:= 0.05, cx0:= Mx0/(1-Mx0),                       Mdb_i:= m_number, ct0:= Mdb_i,                       PDE:= diff(C(x,t),t)=d*diff(C(x,t),x,x),                       IBC:= {C(x,0)=ct0, C(0,t)=cx0, D(C)(L,t)=0};                 return pdsolve(PDE, IBC, numeric);           end proc :   Q:= proc(d)            local pds, i, S:=0:            if   not type(d, numeric)            then return 'procname(_passed)'            end if;            pds:=getSol(d);            for i from 1 to 5 do                pds:-value( t=t_number[i],                            output=listprocedure                          );                  # get the solution at the desired time                eval(C(x,t), %);              # extract the C(x,t) at that time                int(%, 0..2, numeric) / L;    # compute the average of C(x,t) at that time                S += (% - m_number[i])^2;     # accumulate the residuals            end do;            return sqrt(S);       end proc:
 > # # Rouben's calculation # # A few representative values - shows that the # "optimum" value for the parameter 'd' is # somewhere around 0.14 #   Q(0.1), Q(0.2), Q(0.5), Q(1.0);   plot(Q(d), d=0.1 .. 0.3, numpoints=5, view=0..2);  > # # Use the Optimization() package to get the value # of diffusion coefficient which produces the # minimum overall residual. Compute the solution # at this "optimal" diffusion coefficient #   oSol:= Optimization:-Minimize(Q, 0..1);   pdsOpt:= getSol( oSol ): # # Using the optimal diffusion coefficient compute # the average residual at each of the time-values # in the vector t_number. These values *ought* to # match somewhat with the corresponding value in # the vector m_number #   resid:= Vector( 5,                   i-> int                       ( eval                         ( C(x,t),                           pdsOpt:-value                                   ( t=t_number[i],                                     output=listprocedure                                   )                         ),                         0..2,                         numeric                       )/L                  ); # # Get the mean square error in the average C(x,t)-value # at each of the time-values in the vector t_number #   `^`~((resid-m_number),2);   (1)
 > # # Plot C(x,t) for the optimal value of the diffusion # coefficient for each of the time-values in t_number #   cols:= [red, green, blue, brown, black]:   plots:-display( [ seq                     ( pdsOpt:-plot                       ( t=t_number[j],                         x=0..L,                         color=cols[j]                       ),                       j=1..numelems(t_number)                     )                   ]                 ); # # Plot the "general 3-D" solution #   pdsOpt:-plot3d( x=0..L, t=0..t_number[-1]);  >

## Observation...

In this definition, what is s(t)?

Some totally unknown, completely random function perhaps??

Because until it is defined in some way, this system cannot be solved!

While you are thinking about this, you might want to think about the parameter 's0' - now that wouldn't be s(0), would it? Because if you define s(t) using an additional differential equation, then you will need another boundary condition. (Obviously, if it is a straightforward function definition, then you won't, and 's0' will be a 'simple' parameter)

## Any numerical method...

requires that the variable be 'L' be given a value.

The attached uses L=10, althogh it can easily be changed

 > restart; # # Have to give 'L' a value to enable # numeric solution #   L:=10:   PDE := 35139.16048*(diff(w(x, t), x, x, x, x))+98646.00937*(diff(w(x, t), t, t))-2771.636*(diff(w(x, t), x, x)) = 24883.20000:   bcs:= w(0,t)=0, D(w)(0,t)=0, D[1,1](w)(L,t)=0, D[1,1,1](w)(L,t)=0:   ics:= w(x,0)=0, D(w)(x,0)=0:   sol:=pdsolve(PDE, [bcs, ics], numeric): # # Plot w(x,t) for x=L/2 and t=0..10 #   sol:-plot(x=L/2, t=0..10); # # Plot w(x,t) for x=0..L and t=5 #   sol:-plot(t=0.5, x=0..L); # # 3D plot of w(x,t) for x=0..L, and t=0..10) #   sol:-plot3d( x=0..L, t=0..10);   >
 >
 >

## Two possibilities...

You cannot really expect the drawing tool within Maple to have the capabilities of a serious drawing program. Maple is a system for mathemetical calculation - not drawing.

You have two choices

1. Use an external drawing packing (Gimp/Visio/whatever). Export from that package into a format understood by the Maple drawing tool, and then use an "import" command
2. Use the plottoools() package within Maple to produce a suitable drawing, export the result in an appropriate format. Then use an import command in the target worksheet as in (1) above

## Firstly...

upload actual worksheet using the big green up-arrow in the Mapleprimes toolbar.

The code you present in your oriignal post cannot produce the output you present for three reasons

1. the variable 'q' is undefined
2. you have two first order differential equations, which obviously require two initlal conditions in order for a numeric solution to be feasible. No initial conditions are supplied

## Not sure I understand the problem...

What is wrong with a simple function definition, as in the attached; the function may be called from top-level or from within a package. Scoping is pretty much the same as any other function definition

 > restart;   r:= (x,y)-> Record( 'name'=x, 'age'=y):   myModule:= module()                     global r;                     option package;                     export makeRec;                     makeRec:= proc( a, b )                                     return r(a,b)                               end proc;              end module:   with(myModule):
 > s1:=r("JohnDoe1", 26):   s2:=r("JohnDoe2", 28):   s3:=makeRec("JohnDoe3", 30):   s1:-name;   s2:-age;   s3:-age;   (1)
 >
 >

## +ve and -ve...

With R=0.5, the integrand is positive over some of the domain z=-R..R, y=0..R and negative over the rest. My "calibrated eyeball" thnks that, overall, the integral ought to be negative - which it is. See the attached

No idea what your plot command is doing because 'phi' is undefined

 > Integrand := parse(FileTools[Text][ReadFile]("C:/Users/TomLeslie/Desktop/1.txt")): # # Evaluate the integral #   R:=0.5:   evalf(Int(Integrand, [z = -R .. R, y = 0 .. R])); # # Plot the integrand over the region specified # by the limits of integration #   plot3d(Integrand, z=-R..R, y=0..R);  >

## An alternative...

which doesn't use the 'background' option. (and just because I like alternatives)

Rendering is better in the worksheet than it is on this site - honest!

Differential Geometry

 > restart;
 > with(VectorCalculus):with(LinearAlgebra):with(plots):
 > Digits:=10: interface(displayPrecision=10): with(VectorCalculus):BasisFormat(false):

Frenet Frame

 > frenet := proc(p)                local vel, speed, unitT, curv, kappa, unitN, unitB, tors:                vel := diff(p,s):                speed := simplify(Norm(vel)):                unitT := simplify(vel/speed):                curv := simplify(diff(unitT,s)):                kappa := simplify(Norm(curv)):                unitN := simplify(curv/kappa):                unitB := simplify(unitT &x unitN):                tors := diff(unitB,s)/~unitN:               # print("Velocity",vel,"Speed",speed,"T",unitT,"Curvature",curv,"Kappa",kappa,"N",unitN,"B",unitB,"Torsion",tors);                return([vel,curv,unitB]):           end proc:

Testing

 > p :=unapply( <3*arcsinh(s/5), (25+s^2)^(1/2), 4*arcsinh(s/5)>, s): ptan:=unapply(frenet(p(s)),s): display( [ spacecurve            ( p(s),              s = 0 .. 16*Pi,              numpoints = 1000            ),            display            ( [ seq                ( arrow                  ( p(A),                    ptan(A),                    width=0.3,                    length=4                  ),                  A=0..50                )              ],              insequence=true            )          ]        ); >