tomleslie

4923 Reputation

15 Badges

9 years, 187 days

MaplePrimes Activity


These are answers submitted by tomleslie

numbers := [random[empirical[0.5, 0.5]](N)]

uses the 'stats' package, which has been 'deprecated' since Maple version 9.5 in 2004. You cannot access these commands using 'with(Statistics)'. You have to use 'with(stats)' , as shown in the attached - which actually "works".

Note that you really shouldn't be using such deprecated packages without a very very good reason. Since all this command actually does is produce a list containing the integers 1,2 at random - there really doesn't seem to be any sensible reason for continuing to use it!

Also

  1. in the final for loop: when you are just adding a finite list of "known" numbers, it is much more efficient to use add() rather than sum()
  2. the procedure Wpath() is defined but never used - deliberate?


 

restart

with(stats)

Wpath := proc (steps, t) local walk, i, N, ww; N := nops(steps); walk[0] := 0; for i from 0 to N-1 do walk[i+1] := walk[i]+steps[i+1]*sqrt(t/N) end do; ww := seq(plot(walk[i], t*i/N .. t*(i+1)/N), i = 0 .. N-1); plots[display]([ww]) end proc

N := 400

numbers := [random[empirical[.5, .5]](N)]

st1 := map(proc (x) options operator, arrow; 2*x-3 end proc, numbers)

list_of_k := [40, 20, 10, 5, 2, 1]

for j to nops(list_of_k) do k := list_of_k[j]; st[k] := [seq((sum(st1[p], p = i*k-k+1 .. i*k))/sqrt(k), i = 1 .. N/k)] end do

40

 

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

 

20

 

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

 

10

 

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

 

5

 

[-(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), 5^(1/2), (1/5)*5^(1/2), -5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (3/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (3/5)*5^(1/2), -(1/5)*5^(1/2), -5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), -(1/5)*5^(1/2), 5^(1/2), -(1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), -5^(1/2), -(1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), 5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(3/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), 5^(1/2), -5^(1/2), (1/5)*5^(1/2), -(1/5)*5^(1/2), (1/5)*5^(1/2), (1/5)*5^(1/2), (3/5)*5^(1/2), (3/5)*5^(1/2), -(3/5)*5^(1/2), -(3/5)*5^(1/2), 5^(1/2), -(3/5)*5^(1/2), (1/5)*5^(1/2)]

 

2

 

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

 

1

 

[-1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1]

(1)

``


 

Download walk.mw

 

 

Somehow, the "code" in your file Matrix.mw appears to be "Maple Output" which is odd

If I convert everything to 1-D Maple Input, then everything seems to execute correctly - see the attached

  restart;
  interface(rtablesize = 10);

  with(LinearAlgebra);

  A := 8:
  B := 5:
  q := 0.4:
  p := 0.2:
  r := 1 - p - q:
  dimP := A + B + 1:

  P := Matrix(dimP, dimP, [0 $ dimP*dimP]):
  P[1, 1] := 1:
  P[1, 2] := 0:
  P[dimP, dimP] := 1:
  P[dimP, dimP - 1] := 0:
  for i from 2 to dimP - 1 do
      P[i, i - 1] := q;
      P[i, i] := r;
      P[i, i + 1] := p;
  end do:

  p0 := Matrix(dimP, 1, [0 $ dimP]):
  p0[A + 1, 1] := 1:
  pV[0] := p0:
  PT := Transpose(P):

  for n to 200 do
      pV[n] := PT . (pV[n - 1]):
  end do:

  map(x -> evalf(x, 3), Transpose(pV[5]));

_rtable[18446744074332538990]

(1)

 

NULL

Download mat2.mw

what you are getting wrong unless you upload the relevant worksheet using the big green up-arrow in the toolbar.

The attached will produce the plot you want

plot(x^2-9, x=-10..10);

 

 

Download simplePlot.mw

is to get rid of the integral term, by differentiating the pde with respect to 't'.

The drawback with this approach is hte at it increases the degree of the PDE (with respect to 't') so an additional initial condition is needed. I tried a few possibilities for D[2](u)(x,0). These didn't make a huge difference to the overall "form" of the solution, whihc always gets *seriously* huge for t>~1.
You might want to experiment with the following, and in particular with the "extra" initial condition D[2](u)(x,0=1, which I inserted. (Changing the rhs of this initial condition doesn't make a huge difference to the overall form!)

interface(rtablesize=10):
f:=(x^3+t^2*x^2-t^2*x+4*t*x-2*t+1)*exp(x*t)-x+1;
pde:=diff(u(x,t),t)+diff(u(x,t),x,x)+u(x,t)+int(u(x,tau),tau=0..t)=f;
pde:=diff(pde,t);
bounds:= u(0,t)=0, u(1,t)=0, u(x,0)=x*(x-1), D[2](u)(x,0)=1;
psol:=pdsolve(pde, [bounds], numeric);
psol:-plot3d( u(x,t), x=0..1, t=0..0.1);

(t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*exp(t*x)-x+1

 

diff(u(x, t), t)+diff(diff(u(x, t), x), x)+u(x, t)+int(u(x, tau), tau = 0 .. t) = (t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*exp(t*x)-x+1

 

diff(diff(u(x, t), t), t)+diff(diff(diff(u(x, t), t), x), x)+diff(u(x, t), t)+u(x, t) = (2*t*x^2-2*t*x+4*x-2)*exp(t*x)+(t^2*x^2-t^2*x+x^3+4*t*x-2*t+1)*x*exp(t*x)

 

u(0, t) = 0, u(1, t) = 0, u(x, 0) = x*(x-1), (D[2](u))(x, 0) = 1

 

_m906819456

 

 

 

Download pdeSol2.mw

the attached -maybe?

  restart;
  interface(rtablesize=12);
  M2:=Matrix(2, 12, (i,j)->`if`( i=1,
                                 `if`(type(j, odd), U[(j+1)/2],0),
                                 `if`(type(j, even), U[j/2], 0)
                               )
            );
  Matrix( [ [seq( D[1](M2[1,j])(x,y) + 0*M2[2,j],          j =1..12)],
            [seq( 0*(M2[1,j])(x,y)   + D[2](M2[2,j])(x,y), j =1..12)],
            [seq( D[2](M2[1,j])(x,y) + D[1](M2[2,j])(x,y), j =1..12)]
          ]
        ):
  convert~([%], diff);

[10, 10]

 

Matrix(2, 12, {(1, 1) = U[1], (1, 2) = 0, (1, 3) = U[2], (1, 4) = 0, (1, 5) = U[3], (1, 6) = 0, (1, 7) = U[4], (1, 8) = 0, (1, 9) = U[5], (1, 10) = 0, (1, 11) = U[6], (1, 12) = 0, (2, 1) = 0, (2, 2) = U[1], (2, 3) = 0, (2, 4) = U[2], (2, 5) = 0, (2, 6) = U[3], (2, 7) = 0, (2, 8) = U[4], (2, 9) = 0, (2, 10) = U[5], (2, 11) = 0, (2, 12) = U[6]})

 

[Matrix(%id = 18446744074423018070)]

(1)

 


 

Download doDiffs.mw

which appears to give the same answer

restart

N := 2; A := -N; B := N

q := .3; p := .5; sa := .9; sb := .1; r := 1-p-q

dimP := 2*N+1

P := Matrix(dimP, dimP)

P[1, 1] := sa; P[1, 2] := 1-sa; P[dimP, dimP] := sb; P[dimP, dimP-1] := 1-sb

for i from 2 to dimP-1 do P[i, i-1] := q; P[i, i] := r; P[i, i+1] := p end do

P

Matrix(%id = 18446744074421179446)

(1)

# change this part code to the modernversion
 with(linalg)

J := diag(`$`(1, dimP))

array( 1 .. 5, 1 .. 5, [( 5, 5 ) = (1), ( 3, 3 ) = (1), ( 4, 4 ) = (1), ( 2, 2 ) = (1), ( 1, 1 ) = (1)  ] )

(2)

d := matrix(dimP, 1, [`$`(1, dimP)])

array( 1 .. 5, 1 .. 1, [( 3, 1 ) = (1), ( 4, 1 ) = (1), ( 5, 1 ) = (1), ( 1, 1 ) = (1), ( 2, 1 ) = (1)  ] )

(3)

b := matrix(dimP+1, 1, [`$`(0, dimP), 1])

array( 1 .. 6, 1 .. 1, [( 3, 1 ) = (0), ( 4, 1 ) = (0), ( 5, 1 ) = (0), ( 1, 1 ) = (0), ( 6, 1 ) = (1), ( 2, 1 ) = (0)  ] )

(4)

augment(P-J, d); A := transpose(augment(P-J, d))

Matrix(5, 6, {(1, 1) = -.1, (1, 2) = .1, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 1, (2, 1) = .3, (2, 2) = -.8, (2, 3) = .5, (2, 4) = 0, (2, 5) = 0, (2, 6) = 1, (3, 1) = 0, (3, 2) = .3, (3, 3) = -.8, (3, 4) = .5, (3, 5) = 0, (3, 6) = 1, (4, 1) = 0, (4, 2) = 0, (4, 3) = .3, (4, 4) = -.8, (4, 5) = .5, (4, 6) = 1, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = .9, (5, 5) = -.9, (5, 6) = 1})

 

array( 1 .. 6, 1 .. 5, [( 5, 4 ) = (.5), ( 4, 3 ) = (.5), ( 5, 5 ) = (-.9), ( 4, 5 ) = (.9), ( 6, 5 ) = (1), ( 1, 4 ) = (0), ( 3, 3 ) = (-.8), ( 4, 2 ) = (0), ( 3, 1 ) = (0), ( 2, 4 ) = (0), ( 4, 4 ) = (-.8), ( 5, 2 ) = (0), ( 1, 3 ) = (0), ( 1, 5 ) = (0), ( 3, 4 ) = (.3), ( 6, 4 ) = (1), ( 5, 3 ) = (0), ( 2, 2 ) = (-.8), ( 4, 1 ) = (0), ( 1, 2 ) = (.3), ( 5, 1 ) = (0), ( 3, 2 ) = (.5), ( 2, 5 ) = (0), ( 6, 3 ) = (1), ( 1, 1 ) = (-.1), ( 2, 3 ) = (.3), ( 3, 5 ) = (0), ( 6, 2 ) = (1), ( 6, 1 ) = (1), ( 2, 1 ) = (.1)  ] )

(5)

linsolve(A, b)

array( 1 .. 5, 1 .. 1, [( 3, 1 ) = (.1668726823), ( 4, 1 ) = (.2781211372), ( 5, 1 ) = (.1545117429), ( 1, 1 ) = (.3003708282), ( 2, 1 ) = (.1001236094)  ] )

(6)

#
# Unload the linalg package and load the
# LinearAlgebra packages
#
  unwith(linalg):
  with(LinearAlgebra):
#
# Same calculation using LinearAlgebra()
# commands
#
  J:= DiagonalMatrix([1 $ dimP]):
  d:= Matrix(dimP, 1, [1 $ dimP]):
  b:= Matrix(dimP + 1, 1, [0 $ dimP, 1]):
  A:= Transpose(<P-J|d>):
  LinearSolve(A,b);
 

Vector(5, {(1) = .30037082818294225, (2) = .10012360939431399, (3) = .16687268232385657, (4) = .27812113720642756, (5) = .15451174289245975})

(7)

 


 

Download modernize.mw

  1. As defined 'A' is complex
  2. So everything subsequently based on 'A' will also be complex (??)
  3. Doesn't seem to cause an issue for 'H' and 'B' -maybe because imaginary part are very small???
  4. Can't fieldplot 'E', but can fieldplot 'Re(E)'. It is noticeable that for 'E', imaginary parts are "bigger" than real parts

See the attached

  restart;

  with(plots):
  with(LinearAlgebra):
  with(VectorCalculus):
  SetCoordinates('cartesian'[x,y,z]):

  mu_0:=4*Pi*10^(-7): I_0:=2400: f:=2500: omega:=2*Pi*f:
  c:=3*10^8: k:=omega/c: l:=3*10^(-2): epsilon_0:=1/(mu_0*c^2):

  J:=Vector[column]([ 0 ,
                      0 ,
                      I_0/s ]);

Vector(3, {(1) = 0, (2) = 0, (3) = 2400/s})

(1)

# IMPORTANT: R is constant for the calculation of A
#
  R:=sqrt(x^2+y^2+z^2);
#
# Note 'A' will be complex!
#
  A := VectorField(mu_0/(4*Pi)*int(J*exp(-I*k*R)*s/R, p = -l/2 .. l/2));

R := sqrt(x^2+y^2+z^2)

 

Vector(3, {(1) = 0, (2) = 0, (3) = (9/1250000)*exp(-((1/60000)*I)*Pi*sqrt(x^2+y^2+z^2))/sqrt(x^2+y^2+z^2)})

(2)

  B:=Curl(A):
#
# B ought(?) to be complex, because 'A' is - yet
# it plots??? Is fieldplot() perhaps ignoring
# "small" imaginary parts??
#
# No difference between plotting 'B' and 'Re(B)'
#
  B_plot:=fieldplot3d(B,     x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);
  B_plot:=fieldplot3d(Re(B), x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM)

 

 

#
# H is complex, because B is. Maybe fieldplot ia
# still ignoring "small" imaginary parts??
#
# No difference between plotting 'H' and 'Re(H)'
#
  H:=(1/mu_0)*B:
  H_plot:=fieldplot3d(H, x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);
  H_plot:=fieldplot3d(Re(H), x=-1..1,y=-1..1,z=-1..1,fieldstrength=log,arrows=SLIM);

 

 

#
# E is complex, because 'H' is. Now there
# is a significant difference between plotting
# 'E' and 'Re(E)'
#
  E:=1/(omega*epsilon_0*I)*Curl(H):
  E_plot:=fieldplot3d(E, x=1..500,y=1..500,z=1..500,fieldstrength=log,arrows=SLIM);
  E_plot:=fieldplot3d(Re(E), x=1..500,y=1..500,z=1..500,fieldstrength=log,arrows=SLIM);

 

 

#
# Check relative magnitudes of real and imaginary
# parts in B, H, and E
#
  evalf(eval(B, [x=1,y=1,z=1]));
  evalf(eval(H, [x=1,y=1,z=1]));
  evalf(eval(E, [x=1,y=1,z=1]));

Vector(3, {(1) = -0.1385640652e-5+0.5000000000e-18*I, (2) = 0.1385640652e-5-0.5000000000e-18*I, (3) = 0.})

 

Vector(3, {(1) = -1.102657795+0.3978873575e-12*I, (2) = 1.102657795-0.3978873575e-12*I, (3) = 0.})

 

Vector(3, {(1) = -0.8908680955e-6-7939136.102*I, (2) = -0.8908680955e-6-7939136.102*I, (3) = -0.3947841759e-5-0.4353118458e-1*I})

(3)

 


 

Download fields.mw

you check the help at ?solve/details you will find that the solve() command will ignore assumptions which cannot be written as polynomial inequalities. Hence you cannot use 'assumptions' to obtain only real solutions.

For real solutions, the recommendation is to use the RealDomain() package. Alternatively you can 'filter' the roots to select only the 'real' ones.

Both approaches are shown in the attached
 

restart;
with(RealDomain):
fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):
sol := solve(And(diff(g(x),x)=0,diff(g(x),x,x)<>0),x);
evalf(sol);

-1, -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2)

 

-1., -.3051050611, -1.694894939, .497747328, -2.497747328, .283407770, -2.283407770

(1)

restart;
fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):
sol := solve(And(diff(g(x),x)=0,diff(g(x),x,x)<>0),x):
realSol:=seq( `if`( type(evalf(j), realcons), j, NULL), j in sol);
evalf(realSol);

-1, -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 1))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 2))^(1/2), -1+(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2), -1-(1+RootOf(2*_Z^6-3*_Z^5+4*_Z^4+5*_Z^3-14*_Z^2+4, index = 4))^(1/2)

 

-1., .283407770, -2.283407770, .497747328, -2.497747328, -.3051050611, -1.694894939

(2)

 


 

Download realSols.mw

 

just how "fussy" you want to get about the positioning of the labels, adn exactly how they move

The attached is pretty basic, but works

Download rollCube.mw

  with(GraphTheory):

  truss:=Graph( { {1, 2}, {1, 3}, {2, 3},
                  {2, 4}, {3 ,4}, {3, 5},
                  {4, 5}, {4, 6}, {5, 6}
                }
              ):
  SetVertexPositions( truss,
                      [ [0, 1], [0, 0], [1, 1],
                        [1, 0], [2, 1], [2, 0]
                      ]
                    ):

  DrawGraph(truss);

 

 


Download truss.mw

You can adjust vertexColors, edgeColors, etc, etc - if yu really want to

If I replace the plot() command in your code with

seq( evalf(eval(c1,[gamma2=0.02, gamma1=j])), j=0.02..0.1,0.02)

just to test what values might be obtained, then I get

-0.05316538521 - 1.597731966*I,
-0.05642551576 - 1.600665914*I,
-0.05990502879 - 1.603524111*I,
-0.06361988171 - 1.606284242*I,
-0.06758719718 - 1.608919583*I

In other words every value is complex. Is this desired or expected?

If it is "expected", how were you planning to plot it?

Fairly obviously, for this example, if you evaluate 'sol' under the condition 'f(x)=0', you are going to get divide-by-zero errors.

There are a few wayst to get around this problem - one of them given by vv - and another as shown in the attached

restart;

depvars:= (pde::{algebraic, `=`(algebraic), {set,list}(algebraic, `=`(algebraic))})->
   indets(
      indets(convert(pde, diff), specfunc(diff)),
      And(typefunc(name, name), Not(typefunc({mathfunc, identical(diff)})))
   );

is_solution_trivial:= (pde, sol::{`=`, set(`=`),`function`})->
   evalb(eval(`if`(sol::set, sol, {sol}), depvars(pde)=~ 0) = {0=0});

proc (pde::{algebraic, `=`(algebraic), ({list, set})(algebraic, `=`(algebraic))}) options operator, arrow; indets(indets(convert(pde, diff), specfunc(diff)), And(typefunc(name, name), Not(typefunc({mathfunc, identical(diff)})))) end proc

 

proc (pde, sol::{`=`, function, set(`=`)}) options operator, arrow; evalb(eval(`if`(sol::set, sol, {sol}), `~`[:-`=`](depvars(pde), ` $`, 0)) = {0 = 0}) end proc

(1)

pde := diff(w(x,y),x)-( diff(f(x),x)*y^2 -f(x)*g(x)*y+ g(x))*diff(w(x,y),y) = 0;

diff(w(x, y), x)-((diff(f(x), x))*y^2-f(x)*g(x)*y+g(x))*(diff(w(x, y), y)) = 0

(2)

sol:=pdsolve(pde,w(x,y));

w(x, y) = _F1((y*f(x)*(Int((diff(f(x), x))*exp(Int(f(x)*g(x), x))/f(x)^2, x))-f(x)*exp(Int((f(x)^2*g(x)-2*(diff(f(x), x)))/f(x), x))-(Int((diff(f(x), x))*exp(Int(f(x)*g(x), x))/f(x)^2, x)))/(y*f(x)-1))

(3)

is_solution_trivial(pde,sol)

 

Error, (in eval/Int) numeric exception: division by zero

 

#
# Check output of depvars() - returns the names
# of 'dependent' functions which have derivatives
# in the pde (or pde system)
#
# Does not return the names of other 'dependent'
# functions (eg g(x) in this example) which do
# not have derivatives in the pde
#
  depvars(pde);

{f(x), w(x, y)}

(4)

#
# The is_solution_trivial() function
#
#  1. sets up the equations depvars(pde)=~0, So,
#     in this example, one obtains
#  2. evaluates its second argument (in this
#     case 'sol') under the condition depvars(pde)=~0
#     ie {f(x) = 0, w(x, y) = 0}
#  3. Thus (for this example), the whole process is
#     equivalent to
#
#     eval(sol, {f(x) = 0, w(x, y) = 0} )
#
#     Since the term 1/f(x) occurs multiple times
#     in sol, this becomes 1/0 and the division-
#     by-zero error is inevitable.
#
# So the following is also guaranteed to give the
# same error
#
  eval(sol, {f(x) = 0, w(x, y) = 0} );

Error, (in eval/Int) numeric exception: division by zero

 

#
# One could circumvent this issue using something
# like
#
  if   not has(sol, 1/~depvars(pde))
  then is_solution_trivial(pde,sol);
  fi;

 

Download div0.mw

In my Maple 2018

  1. with default setting of Digits, fsolve() cannot find a root, even when it it is given a seriously restricted range over which t search
  2. with default setting of Digits, and no restriction on the possible search range, DirectSearch:-SolveEquations() finds the "correct" root, although the reported residual is "high" - lthough perhaps not "unreasonably" so, given the values which the original expression is capable of generating
  3. setting Digits=30, fsolve() finds the desired root, and DirectSearch:-SolveEquations() finds the same root with negligible residual

See attached

  restart;
  interface(version);
#
# Define parameters
#
  params:= { k = 1.3806e-23,
             T = 300,
             h__bar = 1.05456e-34,
             L__z = 6e-9,
             m__hhw = 3.862216e-31,
             E__hh0 = 3.012136e-21,
             E__hh1 = 1.185628e-20
           }:
#
# Define expression
#
  p:= k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))
      +
      ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z)
      =
      0.3e25;
#
# Plot expression. Note the horizontal plotting range
# and the 'view window' are obtained 'by experiment'.
# But the plot shows that the desired root is somewhere
# around -48
#
  plot( eval(lhs(p)-rhs(p), params), E__fv= -55..-45, numpoints=10000, view=[-55..-45, -1e24..1e24]);

`Standard Worksheet Interface, Maple 2018.2, Windows 7, November 16 2018 Build ID 1362973`

 

k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))+ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z) = 0.3e25

 

 

 

#
# Attempt to obtain the root whilst giving fsolve()
# successively 'tighter' ranges over which to search.
#
# None of these work!
#
  fsolve(eval(p, params));
  fsolve(eval(p, params), E__fv= -infinity..0);
  fsolve(eval(p, params), E__fv= -50..-45);

fsolve(0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.3012136e-20))+0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.1185628e-19)) = 0.3e25, E__fv)

 

fsolve(0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.3012136e-20))+0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.1185628e-19)) = 0.3e25, E__fv, -infinity .. 0)

 

fsolve(0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.3012136e-20))+0.7631009085e25*ln(1+0.2414409194e21*exp(E__fv-0.1185628e-19)) = 0.3e25, E__fv, -50 .. -45)

(1)

#
# Try the equivalent DirectSearch() command with no
# range restriction at all
#
# Looks like this get the "correct" answer, although
# the residual is pretty high - although maybe not
# that high given the values which the expression
# generates! See the y-scale in the above plot
#
  DirectSearch:-SolveEquations(eval(p, params));

[2.334666046828865*10^19, Vector(1, {(1) = -4831838208.}), [`#msub(mi("E"),mi("fv"))` = -48.46001884338331], 69]

(2)

#
# Can reduce the residual in the DirectSearch() by
# increasing the setting of Digits
#
  Digits:=30:
  DirectSearch:-SolveEquations(eval(p, params));

[0., Vector(1, {(1) = 0.}), [`#msub(mi("E"),mi("fv"))` = -48.4600188438469366532693870793], 100]

(3)

#
# With this increased setting of Digits, fsolve()
# (with no range restriction) also works
#
  fsolve(eval(p, params));

-48.4600188438469366532693870793

(4)

 

Download getRoot.mw

 

you must be confusing 'indices' and 'entries'???

Contemplate the attached which takes two lists 'x' and 'y', defines two tables ('a' and 'b').

The first table uses values from 'x' as indices and values from 'y' as entries

The second table uses values from 'y' as indices and values from 'x' as entries

The worksheet also illustrates various ways to show the entries and indices of the tables 'a' and 'b'

#
# Produce a list of more or less random values
#
  x:=[3, 25, 17, 1, 7];
#
# Define another list which has some relation to
# that specified by the OP
#
  y:=[6, 6, 3, 3, 3];
#
# Define a table with indices from the list 'x'
# above, and entries from the list 'y' above
#
  a:=table([seq(x[i]=y[i], i=1..numelems(x))]):
#
# Define another table with indices from the
# list 'y' above, and entries from the list
# 'y' above
#
  b:=table([seq(y[i]=x[i], i=1..numelems(y))]):

[3, 25, 17, 1, 7]

 

[6, 6, 3, 3, 3]

(1)

#
# Check the indices and entries of the table 'a'
#
  indices(a, 'nolist');
  entries(a, 'nolist');

1, 3, 7, 17, 25

 

3, 6, 3, 3, 6

(2)

#
# Check the indices and entries of the table 'b'
#
  indices(b, 'nolist');
  entries(b, 'nolist');

3, 6

 

7, 25

(3)

#
# Show both tables in the form
#
#    indexValue=entryValue
#
  indices(a, pairs);
  indices(b, pairs);

  entries(a, pairs);
  entries(b, pairs);

1 = 3, 3 = 6, 7 = 3, 17 = 3, 25 = 6

 

3 = 7, 6 = 25

 

1 = 3, 3 = 6, 7 = 3, 17 = 3, 25 = 6

 

3 = 7, 6 = 25

(4)

  

 


 

Download tableStuff.mw

You can "see" the behaviour by using the inert form %assuming as  in the attached

restart;
result:=int(x*cos(n*Pi/5*x),x=0..5);
simplify(result) %assuming n::integer and n>0;
value(%);

25*(n*Pi*sin(n*Pi)+cos(n*Pi)-1)/(n^2*Pi^2)

 

%assuming([simplify(25*(n*Pi*sin(n*Pi)+cos(n*Pi)-1)/(n^2*Pi^2))], [false])

 

(25*n*Pi*sin(n*Pi)+25*cos(n*Pi)-25)/(n^2*Pi^2)

(1)

result:=int(x*cos(n*Pi/5*x),x=0..5);
simplify(result) %assuming n::integer, n>0;
value(%);

25*(n*Pi*sin(n*Pi)+cos(n*Pi)-1)/(n^2*Pi^2)

 

%assuming([simplify(25*(n*Pi*sin(n*Pi)+cos(n*Pi)-1)/(n^2*Pi^2))], [n::integer, 0 < n])

 

(-25+25*(-1)^n)/(n^2*Pi^2)

(2)

 

Download assum.mw

 

3 4 5 6 7 8 9 Last Page 5 of 111