tomleslie

4923 Reputation

15 Badges

9 years, 186 days

MaplePrimes Activity


These are answers submitted by tomleslie

- and there are many others!

NB the 3-d plot looks better in a worksheet than it does on this site (which apparently(?) can't deal with the "transparency" option)

  restart;
  with(geom3d):
#
# Cylinder definition
#
  cylR:=3:
  cylH:=2:
#
# Define four points - using symmetry
#
  point(a,  x1,  sqrt(cylR^2-x1^2), 0):
  point(b,  x1, -sqrt(cylR^2-x1^2), 0):
  point(c, -x1, -sqrt(cylR^2-x1^2), cylH):
  point(d, -x1,  sqrt(cylR^2-x1^2), cylH):
#
# determine 'x1' so that length ab is the
# same as length ac
#
  assign(solve([x1>0, distance(b,c)=distance(a,b)], x1));

#
# Get coordinates of a, b, c, d
#
  coordinates~([a, b, c, d])[];
#
# check edge lengths
#
  distance(a, b), distance(b, c), distance(c, d), distance(d, a);
#
# Check diagonals are the same length
#
  distance(a, c), distance(b, d);

[2, 5^(1/2), 0], [2, -5^(1/2), 0], [-2, -5^(1/2), 2], [-2, 5^(1/2), 2]

 

2*5^(1/2), 2*5^(1/2), 2*5^(1/2), 2*5^(1/2)

 

2*10^(1/2), 2*10^(1/2)

(1)

#
# Draw the construction
#
  segment(sab, [a, b]):
  segment(sbc, [b, c]):
  segment(scd, [c, d]):
  segment(sda, [d, a]):
  p1:= draw([sab, sbc, scd, sda], color=black):
  p2:= plottools:-cylinder([0,0,0],3,2, color=red, style=surface, transparency=0.9):
  plots:-display([p1, p2], view=[-5..5, -5..5, -5..5]);

 

 


 

Download sqCyl.mw

when I look at this question, and your previous question posted here

https://www.mapleprimes.com/questions/227697-Help-Plots-Gridlines

I think you might want something like the attached

  restart;
  with(plots):
  with(plottools):
  sol:=solve( [ x^2-y=0, y-(x+5)=0],
              [x, y],
              explicit
            ):
  f:= x-> line([rhs~(x[1]),rhs~(x[2])][], color=red):
  display( [ inequal( [ x^2-y<0,
                        y-(x+5)<0
                      ],
                      x=-4..4,
                      y=0..10,
                      color=white
                    ),
             textplot( [ [ rhs~(sol[1])[],
                           "A",
                           align=right
                         ],
                         [ rhs~(sol[2])[],
                           "B",
                           align=left
                         ]
                       ],
                       font=[times, bold, 18]
                     ),
             seq( f( solve
                     ( [ x^2-y=0, (y+i)-(x+5)=0],
                         [x, y],
                         explicit
                     )
                   ),
                  i=0..5, 0.5
                )
           ],
           size=[600, 600]
         );

 

 

NULL

Download inegplot2.mw

 

 

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
       )

 

 

Download chat.mw

 

 

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[1])[],
                           "A",
                           align=right
                         ],
                         [ rhs~(sol[2])[],
                           "B",
                           align=left
                         ]
                       ],
                       font=[times, bold, 18]
                     )
           ],
           size=[600, 600]
         );

 

 

NULL


 

Download inegplot.mw

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;

x^2*(diff(y(x), x))^2-(1+2*x*y(x))*(diff(y(x), x))+1+y(x)^2 = 0

(1)

DEtools:-odeadvisor(ode)

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

(2)

Vector([dsolve(ode,y(x))]); #now it shows singular solution (first one below)

Vector(3, {(1) = y(x) = (1/4)*(4*x^2-1)/x, (2) = y(x) = x*_C1-sqrt(_C1-1), (3) = y(x) = x*_C1+sqrt(_C1-1)})

(3)

PDEtools:-casesplit(ode)

`casesplit/ans`([(diff(y(x), x))^2 = (2*y(x)*(diff(y(x), x))*x+diff(y(x), x)-y(x)^2-1)/x^2], [2*(diff(y(x), x))*x^2-2*x*y(x)-1 <> 0]), `casesplit/ans`([y(x) = (1/4)*(4*x^2-1)/x], [])

(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)))

Vector(2, {(1) = y(x) = x*(diff(y(x), x))+sqrt(diff(y(x), x)-1), (2) = y(x) = x*(diff(y(x), x))-sqrt(diff(y(x), x)-1)})

(5)

DEtools:-odeadvisor(odes[1]);
DEtools:-odeadvisor(odes[2]);

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

 

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

(6)

dsolve(odes[1],y(x)); #where is singular solution?
dsolve(odes[1],y(x), singsol=all);  #here

y(x) = x*_C1+(_C1-1)^(1/2)

 

y(x) = (1/4)*(4*x^2-1)/x, y(x) = x*_C1+(_C1-1)^(1/2)

(7)

dsolve(odes[2],y(x)); #where is singular solution?
dsolve(odes[2],y(x), singsol=all);  #here

y(x) = x*_C1-(_C1-1)^(1/2)

 

y(x) = (1/4)*(4*x^2-1)/x, y(x) = x*_C1-(_C1-1)^(1/2)

(8)

 

Download singSol.mw

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

 

 

 


 

Download odeSols.mw

 

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):
 

 

[-3/2, -2, 0]

 

[3/2, -2, 0]

 

[3/2, 2, 0]

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

 

Geom3dDetail(["name of the object", tet], ["form of the object", gtetrahedron3d], ["the 4 vertices", [[-3/2, -2, 0], [3/2, -2, 0], [3/2, 2, 0], [-8/3, -3, (1/6)*sqrt(239)]]], ["the 4 faces", [[[-3/2, -2, 0], [3/2, -2, 0], [3/2, 2, 0]], [[-3/2, -2, 0], [3/2, -2, 0], [-8/3, -3, (1/6)*sqrt(239)]], [[-3/2, -2, 0], [3/2, 2, 0], [-8/3, -3, (1/6)*sqrt(239)]], [[3/2, -2, 0], [3/2, 2, 0], [-8/3, -3, (1/6)*sqrt(239)]]]])

(2)

 


Download tetFun2.mw

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/3)*239^(1/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
       );

Matrix(3, 3, {(1, 1) = 2*d[0, 1]^2, (1, 2) = d[0, 1]^2+d[0, 2]^2-d[1, 2]^2, (1, 3) = d[0, 1]^2+d[0, 3]^2-d[1, 3]^2, (2, 1) = d[0, 1]^2+d[0, 2]^2-d[1, 2]^2, (2, 2) = 2*d[0, 2]^2, (2, 3) = d[0, 2]^2+d[0, 3]^2-d[2, 3]^2, (3, 1) = d[0, 1]^2+d[0, 3]^2-d[1, 3]^2, (3, 2) = d[0, 2]^2+d[0, 3]^2-d[2, 3]^2, (3, 3) = 2*d[0, 3]^2})

 

(1/3)*239^(1/2)

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

(1/3)*239^(1/2)

 

[[0, 0, 0], [3, 0, 0], [3, 4, 0], [-7/6, -1, (1/6)*239^(1/2)]]

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

 

[[-3/2, -2, -(99/478)*239^(1/2)], [3/2, -2, -(99/478)*239^(1/2)], [3/2, 2, -(99/478)*239^(1/2)], [-8/3, -3, -(29/717)*239^(1/2)]]

(4)

 

Download tetFun.mw

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?

NULL

Set up functions

 

restart; F := [k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)]; H := alpha*(x[1](t)+x[2](t))
Pars := [k[a1] = 10^(-4), k[a2] = 2*10^(-2), k[a1] = 10^(-4), k[d1] = 10^(-5), k[d2] = 10^(-5), R = 1, k[m] = 10^(-4)]

Start plotting

 

with(plots); FnOut := [op(F), F[1]+F[2]]; plts := NULL; bcs := x[1](0) = 0, x[2](0) = 0, y(0) = 0; for tswitch from 100 by 10 to 200 do Sys := (eval(`~`[`=`]([diff(x[1](t), t), diff(x[2](t), t), diff(y(t), t)], `~`[piecewise](t < tswitch, subs(C[T] = 1, FnOut), subs(C[T] = 0, FnOut))), Pars))[]; p := dsolve([Sys, bcs], numeric); plts := plts, odeplot(p, [t, y(t)], 0 .. 300, size = [800, 800]) end do; display(plts)

[k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t), k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t)+k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)]

 

 

op(F)

k[a1]*C[T]*(R-x[1](t)-x[2](t))-k[d1]*x[1](t), k[a2]*C[T]*(R-x[1](t)-x[2](t))-k[d2]*x[2](t)

(2.1)

NULL


 

Download pwODE.mw

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[1], 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[1](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);

HFloat(2.0125989798253805), HFloat(1.9268252656215228), HFloat(4.344399869795477), HFloat(4.74473396546933)

 

 

#
# 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[2][1] ):
#
# 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);

[.8821785835008983, Vector(1, {(1) = .13886542134771712})]

 

Vector(5, {(1) = 11.499999999999842, (2) = 4.009882865095136, (3) = 1.7332102428767426, (4) = .7664331411209234, (5) = .3558085433772424})

 

Vector[column](%id = 18446744074389251246)

(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]);

 

 

 

Download diffCoeff.mw

 

 

 

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)

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[1](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[2](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);

 

 

 

 

 

 


 

Download solvePDE2.mw

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

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

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;

"JohnDoe1"

 

28

 

30

(1)

 

 


 

Download makeRec.mw

1 2 3 4 5 6 7 Last Page 1 of 110