tomleslie

7350 Reputation

17 Badges

10 years, 345 days

MaplePrimes Activity


These are answers submitted by tomleslie

but if I fix some issues related to variable scope, the attached

pts.mw

will produce the following plot

 

When you capitalise the sum() function, you are making it "inert". The simple change in the code below provides an answer

 restart;
 rho:=1/(2);
 mu:=1/(2);
 T:=1/(4)*(t-x)^(1/(2));
 E[rho,mu](T):=sum((T^(k))/(GAMMA(k*rho+ mu)),k=0..5);
 eq1 := h(t) = (8/3)*(int((t-x)^(1/2)*E[rho, mu](T)*(1/(sqrt(Pi)*x^(1/2))-(3/32)*x^2+(3/16)*exp(-x)-1), t = -1 .. 1));

 

Working correctly the attached

dsol.mw

because you posted a 'picture' of it which is unusable (and I'm not going to retype it.

However you may get something from considering the toy example in the code and attachment below

  restart;
  with(Statistics):
  A := < seq(sin(2*Pi*i/50), i = 1..15) >:
  B := < seq(0.5, i = 1..15) >:
  C := < seq(0.05, i = 1..15) >:
  ErrorPlot(A, xerrors = B, yerrors = C);

which will produce the graph below (although there are many options, bells and whistles which you can include in the ErrorPlot() command

scat.mw

 

 

 

is show below (Actually pretty similar to Kitonum's solution!)

restart;
  L0:=[w,b,g]:
  L:=combinat:-choose([seq(j$4, j in L0)], 4);
  occ:=(el->[seq(s=ListTools:-Occurrences(s,el), s=L0)])~(L);

which produces

 L := [[b, b, b, b], [b, b, b, g], [b, b, b, w], [b, b, g, g], 

  [b, b, g, w], [b, b, w, w], [b, g, g, g], [b, g, g, w], 

   [b, g, w, w], [b, w, w, w], [g, g, g, g], [g, g, g, w], 

   [g, g, w, w], [g, w, w, w], [w, w, w, w]]


     occ := [[w = 0, b = 4, g = 0], [w = 0, b = 3, g = 1], 

       [w = 1, b = 3, g = 0], [w = 0, b = 2, g = 2], 

       [w = 1, b = 2, g = 1], [w = 2, b = 2, g = 0], 

       [w = 0, b = 1, g = 3], [w = 1, b = 1, g = 2], 

       [w = 2, b = 1, g = 1], [w = 3, b = 1, g = 0], 

       [w = 0, b = 0, g = 4], [w = 1, b = 0, g = 3], 

       [w = 2, b = 0, g = 2], [w = 3, b = 0, g = 1], 

       [w = 4, b = 0, g = 0]]



 

Maple will solve your boundary-value problem directly, so converting it to an intial-value problem (the essence of the shooting method) seems rather pointless!

Such conversion from BVPs to IVPs is gnerally only used for ODE systems whic are being solved numerically. (As a general rule numerical IVP methods are "more robust" than numerical BVP methods)

The attached solves your original boundary-value problem, then converts it to to an initial-value problem and solves the latter. Same answer both ways

  restart;
#
# Solve the original system as a boundary-value problem
# ie using the boundary conditions
#
# u(sigma) = 0, u(-sigma) = h
#
  ode := diff((diff(u(y), y) + epsilon*diff(u(y), y)^3)/(1 + epsilon*diff(u(y), y)^2), y) = 0;
  bcs := u(sigma) = 0, u(-sigma) = h;
  sol0 := dsolve([ode, bcs], u(y));
#
# Solve the same system using only one boundary condition,
# namely,
#
# u(-sigma) = h
#
# Note that the solution returned will contain an arbitrary
# constant
#
  bc1:= u(-sigma) = h:
  sol1:= dsolve([ode, bc1], u(y)):
#
# 1. Compute the value of the arbitrary constant which ensures
#    that the second boundary condition above is satisfied
# 2. Differentiate the solution, and evaluate it at the given
#    boundary condition ie u(sigma) = 0, using the value of
#    of the arbitrary constant obtained in (1) above.
#
#  This will generate a further boundary condition of the form
#
#  D(u)(-sigma) = someValue
#
  bc2:= convert(subs(isolate(rhs(eval(sol1, y = sigma)) = 0, _C1), eval(diff(sol1, y), y = -sigma)), D):
#
# Solve the same system with the boundary conditions
#
#  u(-sigma) = 0
#  D(u)(-sigma) = someValue
#
# Note that since both boundary conditions are stated at the
# same value of the independent variable, the ODE system
# is now na "initial value" problem, rather than a "boundary
# value" problem.
#
# Such conversion of BVPs to IVPs, is the essence of the "shooting"
# method
#
# The "new" boundary condtions for the IVP are
#
  bc1;
  bc2;
#
# Solving the same ode with these new boundary conditions
# generates the same answer as before
#
  sol2:= dsolve([ode, bc1, bc2], u(y));
#
# Check
#
  is(sol0=sol2);

shoot.mw

 

 

which really should not be necessary!!

Your problem does look like a *bug* in the UnderlyingGraph() command

  restart;
  with(GraphTheory):
  G1:= Digraph( [a,b,c],
                {[[a,b],2],[[b,c],3],[[c,a],4]}
              );
  Under:= proc(g::Graph)
               uses GraphTheory:
               return Graph( Vertices(g),
                             { seq
                               ( [ convert(j[1], set), j[2] ],
                                 j in Edges(g, weights=true)[]
                               )
                             }
                           ):
          end proc:
  G2:=Under(G1);
  Edges(G1, weights=true);
  Edges(G2, weights=true);

 

the code below

restart;
sol:=solve(0 < M^2 - 4*m^2, M) assuming (0 < M, 0 < m);
plots:-inequal([sol[1][],sol[2][]], M=-10..10, m=-10..10)

whihc produces the plot

 

 

slightly differetn way is shown in the attached

ellip.mw

 

  1. Each of your loops runs from 1 to 6 but the functions U[index], V[index], W[index] are only defined for index =1..2
  2. to perform integrations numerically, it is much more efficient to use evalf(Int(Int(...))), rather than evalf(int(int(...))). Not the capitalisation of the command Int()

The attached generates values for for the 8 integrals which are defined for N=2, reasonably quickly (a few secs on my machine)

someInts.mw

provided that you only want something which "looks pretty"

restart;
for i to 2 do
    for j to 2 do 
        S[i, j] := 2*mu*varepsilon[i, j] - add( (2*mu)/3*varepsilon[r, r]*delta[i, j],
                                                r = 1 .. 3
                                              );
        print
        ( parse
          ( cat("S[", i, j, "]")
          )
          =
          parse
          ( StringTools:-Subs
            ( [ "," = "" ],
              convert
              ( S[i,j],
                string
              )
            )
          )
        ); 
     end do;
end do;

which produces

 

 

 

 

 

is shown  by

  restart:
  R:=10:
  r:=5:
  plot3d( [ ( R+r(u,v)*cos(v) )*sin(u),
            ( R+r(u,v)*cos(v) )*cos(u),
            r*sin(v)
          ],
          u=0..2*Pi,
          v=0..2*Pi,
          style=patchnogrid,
          scaling=constrained,
          coordinateview=[-15..15, -15..15,-5..5],
          axis[3]=[tickmarks=[3, subticks=4]],
          lightmodel=light3,
          viewpoint = circleleft
        );

which produces

 

 

 

are shown in the code below

  restart;
  h:= 0.03;
  C:= 1;
  rho:= 7800;
  omega[m, n] := 222;
  S:=(Matrix([[rho*h*omega[m,n]^2+100,200],[280,rho*h*omega[m,n]^2+444]])).(Vector(1..2,[[A],[B]]))=-C*(Vector(1..2,[555,6665]));
#
# One way
#
  solve([entries(lhs(S)=~rhs(S),nolist)]);
#
# Another way - setting up the problem slightly
# differently
#
  restart;
  with(LinearAlgebra):
  h:= 0.03;
  C:= 1;
  rho:= 7800:
  omega[m, n]:= 222:
  M:=Matrix([[rho*h*omega[m,n]^2+100,200],[280,rho*h*omega[m,n]^2+444]]);
  V:=-C*Vector(1..2,[555,6665]);
  solve(GenerateEquations(M, [A,B], V));

see the attached

sol.mw

 

 

 

 

to the question here

https://www.mapleprimes.com/questions/231383-Hexagon-Inscribed-And-Exinscrit-In-A-Circle,

then the code

  restart;
  with(geometry):
  with(plots):
  _EnvHorizontalName:= x:
  _EnvVerticalName:= y:
  local gamma:
  R:=3.0:
  angs:=Array( 0..5,
               [0, (1/3)*Pi, 3*Pi*(1/4)+1/5, 7*Pi*(1/6)+2/5, 8*Pi*(1/5), 13*Pi*(1/7)]
             ):
#
# Define the origin OO
#
  point(OO, [0, 0]):
#
# Define the circle
#
  circle(C, [OO, R]):
#
# Define points on the circle as specified by the
# entries in the array 'angs'
#
  seq( point( cat(P, j), [R*cos(angs[j]), R*sin(angs[j])]), j=0..5):
#
# Generate the tangents at these points
#
  seq( TangentLine( cat(T, j), cat(P, j), C),
       j=0..5
     ):
#
# Get the points of intersection of the tangents
#
  seq( intersection( cat(Q, j), cat(T, j), cat(T, irem(j+1, 6) ) ),
       j=0..5
     ):
#
# Produce the line segments which form the hexagon
#
  seq( segment( cat(S, j), cat(Q, j), cat(Q, irem(j+1,6))),
       j=0..5
     ):
#
# Produce the lines through adjacent tangent points
#
  seq( line( cat(L, j), [cat(P, j), cat(P, irem(j+1,6))]),
       j=0..5
     ):
#
# Produce the lines through the hexagon vertices
#
  seq( line(cat(V, j), [cat( Q, j), cat(Q, j+3)]), 
       j=0..2
     ):
#
# The OP's points alpha, beta, gamma - not sure what to
# call these!
#
  pts:=[beta, gamma, alpha]:
  seq( intersection( pts[j+1], cat(L, j), cat(L, j+3)), j=0..2):
  line(M, [alpha, gamma]):
#
# Draw the geometric structures
#
  p1:= draw( [ C(color=black),
               seq( cat(S, j)(color=green, thickness=3), j=0..5),
               seq( cat(L, j)(color=magenta, linestyle=spacedash), j=0..5),
               seq( cat(V, j)(color=red), j=0..2),
               M(color=blue, thickness=3)
             ],
             view=[-6..10, -15..6],
             axes=none,
             scaling=constrained,
             size=[800, 800]
           ):
#
# Place text at appropriate locations
#
  p2:= textplot( [ [ coordinates(alpha)[], "alpha", align={above, left}, font=[times, bold,18]],
                   [ coordinates(beta)[],  "beta",  align={above, left}, font=[times, bold, 18]],
                   [ coordinates(gamma)[], "gamma", align={above, left}, font=[times, bold, 18]],
                    seq( [coordinates(cat(P,i))[], cat("P", i),align = [above, left] ], i = 0 .. 5),
                    seq( [coordinates(cat(Q,i))[], cat("Q", i),align = [above, left] ], i = 0 .. 5)
                 ],
                 axes=none
               ):
#
# Display geometry and text
#
  display( p1,p2);

will produce the output

 

I didn't do this in my original response, because I thought it made the plot a bit too "cluttered"

See the attached

hex2.mw

 

MAny of the calculations which you are performing can be performed more siimply by using commands from the geometry package.

The code

  restart;
  with(geometry):
  with(plots):
  _EnvHorizontalName:= x:
  _EnvVerticalName:= y:
  local gamma:
  R:=3.0:
  angs:=Array( 0..5,
               [0, (1/3)*Pi, 3*Pi*(1/4)+1/5, 7*Pi*(1/6)+2/5, 8*Pi*(1/5), 13*Pi*(1/7)]
             ):
#
# Define the origin OO
#
  point(OO, [0, 0]):
#
# Define the circle
#
  circle(C, [OO, R]):
#
# Define points on the circle as specified by the
# entries in the array 'angs'
#
  seq( point( cat(P, j), [R*cos(angs[j]), R*sin(angs[j])]), j=0..5):
#
# Generate the tangents at these points
#
  seq( TangentLine( cat(T, j), cat(P, j), C),
       j=0..5
     ):
#
# Get the points of intersection of the tangents
#
  seq( intersection( cat(Q, j), cat(T, j), cat(T, irem(j+1, 6) ) ),
       j=0..5
     ):
#
# Produce the line segments which form the hexagon
#
  seq( segment( cat(S, j), cat(Q, j), cat(Q, irem(j+1,6))),
       j=0..5
     ):
#
# Produce the lines through adjacent tangent points
#
  seq( line( cat(L, j), [cat(P, j), cat(P, irem(j+1,6))]),
       j=0..5
     ):
#
# Produce the lines through the hexagon vertices
#
  seq( line(cat(V, j), [cat( Q, j), cat(Q, j+3)]), 
       j=0..2
     ):
#
# The OP's points alpha, beta, gamma - not sure what to
# call these!
#
  pts:=[beta, gamma, alpha]:
  seq( intersection( pts[j+1], cat(L, j), cat(L, j+3)), j=0..2):
  line(M, [alpha, gamma]):
#
# Draw the geometric structures
#
  p1:= draw( [ C(color=black),
               seq( cat(S, j)(color=green, thickness=3), j=0..5),
               seq( cat(L, j)(color=magenta, linestyle=spacedash), j=0..5),
               seq( cat(V, j)(color=red), j=0..2),
               M(color=blue, thickness=3)
             ],
             view=[-6..10, -15..6],
             axes=none,
             scaling=constrained,
             size=[800, 800]
           ):
#
# Place text at appropriate locations
#
  p2:= textplot( [ [ coordinates(alpha)[], "alpha", align={above, left}, font=[times, bold,18]],
                   [ coordinates(beta)[],  "beta",  align={above, left}, font=[times, bold, 18]],
                   [ coordinates(gamma)[], "gamma", align={above, left}, font=[times, bold, 18]]
                 ],
                 axes=none
               ):
#
# Display geometry and text
#
  display( p1,p2);

will produce the same plot as your worksheet.

See the attached worksheet

hex.mw

 

 

 

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