tomleslie

13876 Reputation

20 Badges

15 years, 166 days

MaplePrimes Activity


These are answers submitted by tomleslie

the code

  restart;
#
# Function whihc does the "work"
#
  s:= p-> indets(p)[]=rhs~( `union`
                            ( isolve
                              ( numer
                                ( normal
                                  ( lhs(p)
                                  )
                                )
                              )
                            )
                          ):
#
# The equations
#
  eqs:= [ 1/(x^2+8*x+9)  + 5/(x^2+8*x+5)  - 8/(x^2+8*x+3)  = 0,
          3/(x^2+5*x+10) + 7/(x^2+5*x+2)  - 3/(x^2+5*x+3)  = 0,
          5/(x^2+6*x-7)  + 4/(x^2+6*x+2)  - 7/(x^2+6*x+1)  = 0,
          2/(x^2+11*x+6) + 6/(x^2+11*x-2) - 7/(x^2+11*x+3) = 0
        ]:
#
# Get the integer solutions
#
  ans:=s~(eqs);

will return

ans := [x = {-7, -5, -3, -1}, x = {-4, -3, -2, -1}, x = {-5, -4, -2, -1}, x = {-10, -9, -2, -1}]

 

there may be more than one way to do it!

I used the system as defined in your original post. Your subsequent "answer" has redefined the system somewhat, and I haven't got th time right now to figure out to figure out whether or not the two systems are "equivalent".

You should be aware that the although the figure I generate in wallball2.mw  and that in your answer  are very similar, there is a slight timeshift between the two. I haven't got the time/energy to track down this discrepancy

However the "events" specification I used is substantially different - again I haven't got the time/energy to work out whether the two "events" specifications are in fact equivalent

  restart:
  with(plots):
  sys:= { diff(v(t), t)=s(t)*b(t), 
          diff(x(t), t)=v(t)*b(t),
          s(t)=10*cos(t)-x(t),
          v(0)=0,
          x(0)=0,
          b(0)=1
        }:
  evs:= [ [ [ increasing(x(t)-2), And(s(t)>0)], b(t)=0],
          [ [ increasing(x(t)-2), And(s(t)<0)], v(t)=-v(t)],
          [ increasing(s(t)), b(t)=1],
          [ decreasing(s(t)), b(t)=-1],
          [ [ decreasing(x(t)), And(s(t)<0)], b(t)=0],
          [ [ decreasing(x(t)), And(s(t)>0)], v(t)=-v(t)]
        ] :
  esol:= dsolve( sys,
                 numeric,
                 events=evs,
                 range=0..50,
                 discrete_variables=[b(t)]
               ):
  odeplot( esol,
           [ [t, x(t)],
             [t, b(t)]
           ],
           t=0..50,
           color=[red, blue],
           refine=5,
           gridlines=true, 
           size=[1000, 400]
         )


  animate( pointplot,
           [ ['rhs(esol(t)[4])',0],
             symbol=solidcircle,
             color=red,
             symbolsize=20
           ],
           t=0..50,
           frames=200,
           axes=boxed
         );

which produces the plot

 

and the animation

 

 

Still not doing your own homework I see!!

Well, as a hint,  a couple of Maple one-liners will tell you that there for the first 1000 odd natrural numbers, 302 are prime and 32 are triangular . The `union`() command will tell you that 333 are triangular OR prime, because one number is both triangular AND prime. The `intersect()` command will tell you which number is common to both sets - spoiler alert, it is 3

a case of reversing the velocity when the ball hits either wall?

So the event "specification" would be events=[ [x(t)-2, v(t)=-v(t)], [x(t), v(t)=-v(t)], which reverses the velocity at x(t)=0 and x(t)=2. So the code

  restart;
  with(plots):
  s := t -> 10*cos(t)-x(t):
  sys := { diff(v(t), t)=s(t), 
           diff(x(t), t)=v(t), 
           v(0)=0, 
           x(0)=0
         }:
  sol:=dsolve(sys, numeric, events=[ [x(t)-2, v(t)=-v(t)], [x(t), v(t)=-v(t)]]);
  odeplot( sol, [ [t, x(t)], [t, v(t)]], t=0..5, numpoints=1000, color=[ red, blue]);

will produce the plot

Note that whenever position (in red) hits 0 or 2, the velocity (in blue) reverses sign

wallball.mw

 

Your original  pde is fourth order in 'x' so you need four boundary conditions on 'x' - which you have

Your original pde is second order in 't', so you need two initial conditions on 't' - which you don't have

In the attached, pdeProb.mw    I have "made up" a second initial condition, and the code shown executes correctly

restart;
In := 5.75*10^(-12):            
M := 1000:
E := 10^12:
V := 100:
m := 2300:
K := 10^9:
A := E*In:
B := M*V^2:
C := 2*M*V^2:
F := M + m:
G := K:
tmax := 0.05;
xmin := 0;
L := 1;

pde := A*diff(w(x, t), x $ 4) + B*diff(w(x, t), x $ 2) + C*diff(w(x, t), x, t) + F*diff(w(x, t), t $ 2) + G*w(x, t);
bc := w(0, t) = 0, w(L, t) = 0, D[1, 1](w)(0, t) = 0, D[1, 1](w)(L, t) = 0:

ic := w(x, 0) = 0, D[2](w)(x,0)=1:

pdsA := pdsolve(pde, eval({bc, ic}, L = 1), numeric);

pdsA:-plot3d( w(x,t), x=0..1, t=0..1);

Note that a non-zero choice for the IC which I added is more or less mandatory in order to ensure that a non-zero solution is obtained. Your PDE and all other BCs/ICs is satisfied by w(x,t)=0

 

 

 

 

 

 

From the Maple help for ?pdsolve/numeric

Calling Sequence

pdsolve(PDEsys, conditions, numeric, vars, options)

Parameters

PDEsys - single or set or list of time-dependent partial differential equations in two independent variables

Your system contains functions such as Tg(t, z, r), whihc has three independent variables

I thiink  you now hold the record for the most syntax errors in any worksheet posted here.

So many that it is impossible to work out what your intention is, so in the attached I just maded lots and lots of wild guesses

I suggest that as a minimum you learn

  1. The difference between an equation (eq a=b) and an assignment (eg a:=b). And when you are using an assignment, don't put a space between ':' and '=' when entering ':='  !!!!
  2. The diferent meanings/uses of round and square brackets - '()' and '[]' are not interchangeable

The attached is syntactically correct, but I doubt that it will do what you desire (or indeed anythng useful

egev.mw

when I wish that I could receive $0.01 for every time I have to fix a problem caused by someone using 2D-input in Maple, when they have no understanding of how to use 2D-input.

*********************
RANT ON

The most common problem with users of 2-D input is that they think "whitespace" is just that - ie "whitespace", which does nothing. In fact, with 2D-input in Maple, that is *sometimes* true and *sometimes* not.

If you do not understand the subtleties of 2D input mode, don't use it!!!

RANT OFF
**********************

The attached worksheet

simPlt.mw

produces the plots

and

 

 

 

 

 

 

 

 

 

There are a *lot* of defined scientific constants. If you assigned all of these constants just by "loading the package" you might overwrite some other assignments which exist in your worksheet - that's a very high risk strategy!

With Maple's calling method, you can access any constant you want,  and assign it to any name you want (which might be convenient).

Note also it is probably a good idea to use uneval quotes in a call to (say) the Constant() command.

Consider the code

  restart;
  with(ScientificConstants):
  c:=GetValue(Constant('c'));
  myValueOfc:=GetValue(Constant('c'));
  c:=GetValue(Constant('c'))*GetUnit(Constant('c'));

 

This worksheet

mkGraph.mw

containing the code

  restart;
  with(StringTools):
  with(GraphTheory):
#
# Read data line by line: each line being read
# as a string. (Strip the ":" character,
# assuming one exists.
#
  fname:="C:/Users/TomLeslie/Desktop/infile.txt":
  k:=0:
  while true do
        k:=k+1:
        line:= readline(fname);
        if   type(line, numeric)
        then break
        else l[k]:= Remove( ":", line):
        fi:
  od:
#
# Create the adjacency matrix
#
  M:=Matrix( k-1, k-1, fill=0):
  for i from 1 by 1 to k-1 do
      for j in parse~(Split(l[i]))[2..-1] do
          M[i,j]:=1;
      od:
  od:
#
# Create/draw the graph from the adjacency matrix
#
  G:=Graph(M):
  DrawGraph(G);

will read a text file in the format you supplied (I used the attached -> infile.txt)

You willl have to convert the variable 'fname' in the above to a filePath which is relevant for your machine.

It converts the input to an adjacency matrix, from which it creates a graph (and draws it - see below)

 

The code should(?) work for any supplied file in more-or-less the same format

 

 

 

See the code in the attached, whihc corrects the error


Error, (in PW) cannot determine if this expression is true or false: (1/4)*10^(1/2) < 0.1000000000e-3

However, I keep getting the error you trap , ie "Zero eigenvalue found" - and I have tried varying the input to the procedure.

By the way you don't need a the procedure maxind(). This functionality is provided by the max[index]() command

PWWcorr.mw

since a quintic has six coefficients and you only have five conditions there are an infinite number of solutions. One cn generate these parametrically.

The code below

  restart:
#
# Define system
#
  f:= x->a__5*x^5 + a__4*x^4 + a__3*x^3 + a__2*x^2 + a__1*x + a__0:
  sys:= [ f(1)=2, f(2)=3, f(3)=4, f(4)=6, f(5)=10 ]:
#
# Solve system, regarding a__5 as a parameter. This choice of
# parameter is arbitrary
#
  par:= a__5:
  sol:= solve
        ( sys,
          `union`(indets~(sys, name)[]) minus {par}
        );
#
# Produce a plot for nVals+1 values of the free parameter
# par
#
  nVals:=20:
  plot
  ( [ seq
      ( eval
        ( eval
          ( f(x), sol),
            par=j
        ),
        j=-nVals..nVals
      )
    ],
    x=0..5,
    color=ColorTools:-HueSpread("Orange", nVals+1, 1/(nVals+1) ),
    view=[0..5, -100..100]
  );

produces the family of plots

polyPlot.mw

 

depending on what exactly you want, which is not speicied clearly

The code

Ennep := <u - u^3/3 + u*v^2, v - v^3/3 + v*u^2, u^2 - v^2>:
plot3d(Ennep, u = -2 .. 2, v = -2 .. 2, shading = zhue, lightmodel = light1, orientation = [89, 54]);
plots:-contourplot(convert(Ennep, list), u=-2..2, v=-2..2);
plots:-contourplot3d(convert(Ennep, list), u=-2..2, v=-2..2)

will produce the plots

 

All of these can be produced with various options - you just have to specify which kind of plot you require

I came up with the attached

fourCorr.mw

which produces the "stem" plot

and the "reconstructed" plot

 

 

The original command

Sum(Int(F(n.x), x=-3.5..1.5), n=1..9);

As the Maple help states

To add a finite sequence of values, rather than compute a formula, use the add command.

You are not interested in computing some kind of closed-form "formula" - so don't use sum/Sum. Hence the the firt stage in optimisng the process is

add(Int(F(n.x), x=-3.5..1.5), n=1..9);

The next stage depends very much on whether int(F(n.x), x=-3.5..1.5) has a symbolic solution or not. Since you repeatedly refuse to post the integrand here,  no-one has any idea. If the integral can be solved "analytically", then the optimumal approach is (probably?)

myInt:=unapply( (int(F(n.x), x=-3.5..1.5), n);
add(Int(myInt(n), n=1..9);

If the integral cannot be performed analytically, then the "next best approach" is probably

add(evalf( Int(F(n.x), x=-3.5..1.5)), n=1..9);

This will try a variety of numerical integration methods whihc may or may not be appropriate. Dependng on what you know about the integrand, it may be useful to specify a numerical integration method, something like

add(evalf( Int(F(n.x), x=-3.5..1.5, method=_d01jac)), n=1..9);

probably isn't a bad starting point.

Using the last of the above commands your 'bonus question" would become something like

G:=t-> local n;
           3*add(evalf( Int(f(n,x,t), x=-3.5..1.5, method=_d01ajc)), n=1..4)
           +
           8*3*add(evalf( Int(f(n,x,t), x=-3.5..1.5, method=_d01ajc)), n=5..9):

So for a "toy" example, the code

restart;
f:=(a, b, c)-> a*sin(b*c)+b*sin(a*c):
G:=t-> local n;
       3*add(evalf( Int(f(n,x,t), x=-3.5..1.5, method=_d01ajc)), n=1..4)
       +
       8*3*add(evalf( Int(f(n,x,t), x=-3.5..1.5, method=_d01ajc)), n=5..9):
G(1);

will return -991.6997991

 

 

 

 

 

 

First 51 52 53 54 55 56 57 Last Page 53 of 207