tomleslie

4465 Reputation

10 Badges

9 years, 93 days

MaplePrimes Activity


These are answers submitted by tomleslie

the attached.

Note that (as Carl has pointed out), the firts two operations yield the same result

  EQ := V^3-R*T*V^2/P-(B^2+P*B*R*T*sqrt(T)/(P-A))*V-A*B/(P*sqrt(T)) = 0;
  diff(EQ,T,V);
  diff(EQ,V,T);
  diff(EQ,V,V);

V^3-R*T*V^2/P-(B^2+P*B*R*T^(3/2)/(P-A))*V-A*B/(P*T^(1/2)) = 0

 

-2*R*V/P-(3/2)*P*B*R*T^(1/2)/(P-A) = 0

 

-2*R*V/P-(3/2)*P*B*R*T^(1/2)/(P-A) = 0

 

6*V-2*R*T/P = 0

(1)

 

 


 

Download doDiffs.mw

Command has been "repositioned".

It is now under the "Evaluate" tab

No idea why.

You have defined T[1], T[2] etc as functions. If you wish to evaluate these for some argument, then you need to use T[1](argumentValue). For example

T[1](3);
T[2](7);
or even
T[1](z)

Please stop posting "pictures" of code. Post editable/executable code by uploading your worksheet using the big green up-arrow in the Mapleprimes toolbar

There are several syntax errors, which essentially prevent your worksheet from executing. I have fixed these in the attached.

However (I am nearly certain that) there are further logical errors, which means that this code will not produce what you wish. You need to consider

  1. Do you expect the local variable 'n' in the procedure and the variable 'n' used at the top-level to be the same entity - because they are not! If you wish to evaluate the variable sigma_xx at the value n=loop_index, then the simplest way is probably to change the loop index to 'j' and write eval( sigma_xx, n=j). I have made this change in the second execution group of the attached.
  2. You have two for loops - both of these assign something to Fx, Should one of these be Fy? I have made this change in the second execution group of the attached
  3. Each time either of these for loops is executed, it will overwrite the value obtained on the previous iteration. I have adjusted the code so that  each for loop generates a sequence of values for Fx and Fy
  4. the return statement needs to be modified to return both sequences.

With these changes (some/all of which may not be what you actually want, becuase I'm guessing!), the code executes, and produces values for some of the integrals, although some are 'undefined'. I don't propose to work out why the integration sometimes return undefined values until/unless you confirm the validity of the changes I have made

restart;
F0:=proc(sigma__xx,N)
local  x,y,Fx,Fy: #removed superfluous comma
assume (w,real,w>0):assume (h,real,h>0):

for n from 0 by 2 to N do Fx:=integrate(sigma__xx*cos(w*Zeta*h),Zeta=0..infinity):
end do;
for n from 1 by 2 to N do Fx:=integrate(sigma__xx*sin(w*Zeta*h),Zeta=0..infinity):
end do;
return [Fx]:
end proc;

sigma__xx := -(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n); #added closing parenthesis (may be in wrong place!)

F0(sigma__xx,4); # 'for' loops in proc need numeric final value


 

Warning, `n` is implicitly declared local to procedure `F0`

 

proc (sigma__xx, N) local x, y, Fx, Fy, n; assume(w, real, 0 < w); assume(h, real, 0 < h); for n from 0 by 2 to N do Fx := integrate(sigma__xx*cos(w*Zeta*h), Zeta = 0 .. infinity) end do; for n by 2 to N do Fx := integrate(sigma__xx*sin(w*Zeta*h), Zeta = 0 .. infinity) end do; return [Fx] end proc

 

-(((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n

 

[int((-(((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n-2*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n^2-3*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n)*sin(w*Zeta*h), Zeta = 0 .. infinity)]

(1)

restart;
F0:= proc(sigma__xx,N)
          local  j, x, y, Fx:=NULL, Fy:=NULL:
          assume (w,real,w>0):
          assume (h,real,h>0):
          for j from 0 by 2 to N do
              Fx:= Fx,
                   integrate
                   ( eval
                     ( sigma__xx,
                       n=j
                     )*cos(w*Zeta*h),
                     Zeta = 0..infinity
                   ):
          end do;
          for j from 1 by 2 to N do
              Fy:= Fy,
                   integrate
                   ( eval
                     ( sigma__xx,
                       n=j
                     )*sin(w*Zeta*h),
                     Zeta = 0..infinity
                   ):
          end do;
          return [Fx], [Fy]:
     end proc:

sigma__xx := -(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n):

F0(sigma__xx,5);

[undefined, undefined, (1/16)*exp(-w*h)*Pi*(3*h^4*w^4+32*h^4*w^2-18*h^3*w^3-32*h^3*w+9*h^2*w^2+9*h*w+9)/h^4], [undefined, -(5/24)*w*exp(-w*h)*Pi*(2*h^2*w^2-9*h*w+3)/h^2, (1/240)*w*exp(-w*h)*Pi*(14*h^4*w^4+400*h^4*w^2-105*h^3*w^3-600*h^3*w+70*h^2*w^2+105*h*w+105)/h^4]

(2)

 

 


 

Download doInt.mw

is that you want to compare the results of solving the ODE using Maple's numeric solver with those obtained from your numeric process. And compare both of these with Maple's exact solution.

I have modified your code to do this in the attached. The Maple numeric solver (on default accuracy settings) is slightly more accurate than your numeric process - but not by much!

restart;
#Differential Equation: 
#y''(t)+a(t)y'+b(t)y=f(t)
#IC's:
#y(0)=alpha,y'(0)=beta

k:=2:  
M:=4:  
N:=2^(k-1)*M:
#  Let's define Legendre Polynomials
P[0]:= t -> 1:
P[1]:= t -> t:
for m from 2 to M-1 do
   P[m]:= unapply(expand((2*m-1)*t*P[m-1](t)/m - (m-1)*P[m-2](t)/m), t)
end do:

 
L:=proc(n,m) local a,b;
  a := (2*n-2)/2^k;
  b := 2*n/2^k;
  return piecewise(a <= t and t <= b, P[m](2^k*t-2*n+1)*sqrt(((m+1)/2)*2^k))
end proc:

Ld := Vector(N):
LL := Matrix(N, N):
T := Vector(N):
for j from 1 to N do T[j] := (j-1/2)/N end do:
for n from 1 to 2^(k-1) do
  for m from 0 to M-1 do
i := n+2^(k-1)*m;
Ld[i] := L(n,m)
     end do
end do:
for i from 1 to N do
  for j from 1 to N do
    LL[i, j] := eval(Ld[i], t = T[j])
  end do
end do:

 
pn:=proc(i,n)
  if n=1 then
    return int(Ld[i],t)  
  else
    return int(pn(i,n-1),t)  
  end if
end proc:
##
RHS:= t->f(t)-a(t)*beta-b(t)*(t*beta+alpha):
f:=t->0:
a:=t->1/t:
b:=t->1:
alpha:=1:
beta:=0:
 

R:=Vector(N):
TMP:=Matrix(N,N):
C:=Matrix(N,N):
for i from 1 to N do
  R[i] := evalf(RHS(T[i]));
  tmp := a(t)*pn(i,1)+b(t)*pn(i,2);
  for j from 1 to N do
    TMP[i,j]:=eval(tmp, t = T[j])
  end do
end do:
use LinearAlgebra in
  C := Transpose(LinearSolve(Transpose(LL+TMP), R))
end use:
sol := beta*t+alpha+add(C[m0]*pn(m0,2), m0 = 1 .. N):
y:=unapply(sol,t):
ode:=t*diff(Y(t),t,t)+diff(Y(t),t)+t*Y(t)=0;  
 mExact:=unapply(rhs(dsolve({ode,Y(0)=1,D(Y)(0)=0})),t);
#
# Solve system using Maple numeric methods
# with all options on defaults
#
  mNumer:=eval( Y(t),
                dsolve
                ( {ode,Y(0)=1,D(Y)(0)=0},
                  numeric,
                  output=listprocedure
                )
              );

t*(diff(diff(Y(t), t), t))+diff(Y(t), t)+t*Y(t) = 0

 

proc (t) options operator, arrow; BesselJ(0, t) end proc

 

proc (t) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](t) else _xout := evalf(t) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.10095317511683091e-1, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = -.5}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, t, Y, YP) option `[Y[1] = Y(t), Y[2] = diff(Y(t),t)]`; if t = 0 then if abs(Y[2]) <= 0. then YP[1] := 0; YP[2] := -(1/2)*Y[1] else error "system with provided initial conditions is singular" end if else YP[1] := Y[2]; YP[2] := (-t*Y[1]-Y[2])/t end if; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, t, Y, YP) option `[Y[1] = Y(t), Y[2] = diff(Y(t),t)]`; if t = 0 then if abs(Y[2]) <= 0. then YP[1] := 0; YP[2] := -(1/2)*Y[1] else error "system with provided initial conditions is singular" end if else YP[1] := Y[2]; YP[2] := (-t*Y[1]-Y[2])/t end if; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..2, {(1) = 0., (2) = 1.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then error "initial conditions cannot be changed for systems with removable singularities"; _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(1..3, {(1) = 18446744074581343590, (2) = 18446744074581343854, (3) = 18446744074581344030}), (3) = [t, Y(t), diff(Y(t), t)], (4) = []}); _solnproc := _dat[1]; _pars := map(rhs, _dat[4]); if not type(_xout, 'numeric') then if member(t, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(t, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(t, ["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(t, 'string')); if type(_res, 'list') then return _res[2] else return NULL end if elif member(t, ["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(t, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(t), 'string') = rhs(t); if lhs(_xout) = "initial" then if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 2, rhs(_xout)]) end if elif not type(rhs(_xout), 'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 2, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[2] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(t), 'string') = rhs(t)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(t) else _ndsol := 1; _ndsol := `tools/gensym`("Y(t)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"), _Inert_EXPSEQ(ToInert(_ndsol), _Inert_VERBATIM(pointto(_dat[2][2])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol), _Inert_EXPSEQ(ToInert(t)))) end if end if; try _res := _solnproc(_xout); _res[2] catch: error  end try end proc

(1)

#
# Print all results to screen.
#
# Exact Maple solution
# Numerical Maple solution
# OP's Numerical solution
#
# Errors from Maple's numerical solution are slightly
# less than those obtained from OP's numerical method
#

  printf( "\n\n\n      Time    Maple Exact      "
          "  Maple Numeric       Maple Error    "
          "  OP Numeric          OP ERROR \n\n");
  seq( printf( " %8a %18.15f %18.15f %18.15f %18.15f %18.15f \n",
                 k,
                 evalf[12](mExact(k)),
                 mNumer(k),
                 mNumer(k)-mExact(k),
                 y(k),
                 evalf[12](mExact(k))-y(k) ),
       k=0..1, 1/10
     );




      Time    Maple Exact        Maple Numeric       Maple Error      OP Numeric          OP ERROR

        0  1.000000000000000  1.000000000000000  0.000000000000000  1.000000000000000  0.000000000000000
     1/10  0.997501562066000  0.997501562066904 -0.000000000033096  0.997501576700000 -0.000000014599130
      1/5  0.990024972240000  0.990024970551038 -0.000000001648962  0.990024978300000 -0.000000006101050
     3/10  0.977626246538000  0.977626243360794 -0.000000003139206  0.977626250400000 -0.000000003864050
      2/5  0.960398226660000  0.960398229592673  0.000000002892673  0.960398229600000 -0.000000002925130
      1/2  0.938469807241000  0.938469808879848  0.000000001679848  0.938469801300000  0.000000005880000
      3/5  0.912004863497000  0.912004845584140 -0.000000017915860  0.912004921400000 -0.000000057903760
     7/10  0.881200888607000  0.881200895274771  0.000000006674771  0.881200993000000 -0.000000104400960
      4/5  0.846287352750000  0.846287315917441 -0.000000036882559  0.846287497900000 -0.000000145100960
     9/10  0.807523798123000  0.807523793955245 -0.000000004144755  0.807523978900000 -0.000000180833760
        1  0.765197686558000  0.765197623203748 -0.000000063396252  0.765197891200000 -0.000000204600000

 

#
# Plot errors produced by Maple's numeric method
# and those produced by OP's numeric method
#
   OPErrors:=[seq([k,mExact(k)-y(k)], k=0..1, 1/10)]:
   MapleErrors:= [seq([k,mExact(k)-mNumer(k)], k=0..1, 1/10)]:
   plot( [OPErrors, MapleErrors], color =[red, blue], legend=["OP", "Maple"], axes=boxed);

 


Exact:=[seq([k,mExact(k)], k=0..1, 1/10)]:
Legendre:=[seq([k,y(k)], k=0..1, 1/10)]:
Error:=[seq([k,mExact(k)-y(k)], k=0..1, 1/10)]:
plot([Exact, Legendre, Error], color=[green,blue,red], legend=["Exact", "Legendre", "Error"], axes=box)

 

 


 

Download ODEnumer.mw

Use the menu setting

Tools->Options->Precision

and uncheck the box labelled "Limit expression length to"

I recomment that you use this option with care.

I would display the output here, but this site won't render it properly.

In a Maple worksheet I get a reasonably compact display of 10X6 graphs which fits on a single screen
 

in the following

  1. I have removed a lot of "stuff" which you don't really need
  2. In order to produce numeric values ( essential for plotting), I have assumed that 'm' is the electronic mass and 'hbar' is Planck's constant divided by 2*PI
  3. The obtained solution is complex (ie it contains I=sqrt(-1)). In ordert o guarantee evaluation to numeric vaues, the attached shows plots of abs( f(x,y,3) ), Re( f(x, y, 3) ) and Im( f(x, y, 3) )

Schrdiff(o(t), t, t)dinger PDE on (x,y,t) with initial and boundary conditions. Zero potential: problem number 169

 

Here is the problem 169 specification and solution from the Nasser list:

  restart;
  interface(rtablesize=10):
  interface(version);
  Physics:-Version();

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1.1)

  pde:= I*diff(f(x,y,t),t)=-hbar^2*(diff(f(x,y,t),x$2)+diff(f(x,y,t),y$2))/(2*m);
  ic:= f(x,y,0)=sqrt(2)*(sin(2*Pi*x)*sin(Pi*y)+sin(Pi*x)*sin(2*Pi*y));
  bc:= f(0,y,t)=0,
       f(1,y,t)=0,
       f(x,1,t)=0,
       f(x,0,t)=0;
  sol:= pdsolve([pde,ic,bc]);

I*(diff(f(x, y, t), t)) = -(1/2)*hbar^2*(diff(diff(f(x, y, t), x), x)+diff(diff(f(x, y, t), y), y))/m

 

f(x, y, 0) = 2^(1/2)*(sin(2*Pi*x)*sin(Pi*y)+sin(Pi*x)*sin(2*Pi*y))

 

f(0, y, t) = 0, f(1, y, t) = 0, f(x, 1, t) = 0, f(x, 0, t) = 0

 

f(x, y, t) = 2^(1/2)*sin(Pi*x)*exp(-((5/2)*I)*hbar^2*t*Pi^2/m)*(2*sin(Pi*y)*cos(Pi*x)+sin(2*Pi*y))

(1.2)

#
# Note: above solution is complex although
# the imaginary part is very small (see
# third plot)
#
# Subsequent plots give abs(), Re() and Im()
#
  plot3d
  ( eval
    ( abs(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "abs( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );
 plot3d
  ( eval
    ( Re(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "Re( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );
 plot3d
  ( eval
    ( Im(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "Im( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );

 

 

 

 

Download schPDE.mw

is very similar to that for your previous question

Check the attached

  restart;
  interface(rtablesize=10):
  interface(version);
  Physics:-Version();

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1)

  N:=8;
  X:=Matrix( N+1,
             N+1,
             (i, j)-> x[i-1]^(j-1)
           );

N := 8

 

Matrix(%id = 18446744074377062150)

(2)

 


 

Download doMat2.mw

 - I think??

Anyhow, check out the attached

  restart;
  interface(rtablesize=10):
  interface(version);
  Physics:-Version();
  

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1)

  N:= 8:
  A:= Matrix
      ( N+1,
        N+1,
        (i,j)-> `if`
                 ( j>=i,
                   binomial(j-1, i-1)*(-tau)^(j-i),
                   0
                 )
      );

Matrix(9, 9, {(1, 1) = 1, (1, 2) = -tau, (1, 3) = tau^2, (1, 4) = -tau^3, (1, 5) = tau^4, (1, 6) = -tau^5, (1, 7) = tau^6, (1, 8) = -tau^7, (1, 9) = tau^8, (2, 1) = 0, (2, 2) = 1, (2, 3) = -2*tau, (2, 4) = 3*tau^2, (2, 5) = -4*tau^3, (2, 6) = 5*tau^4, (2, 7) = -6*tau^5, (2, 8) = 7*tau^6, (2, 9) = -8*tau^7, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = -3*tau, (3, 5) = 6*tau^2, (3, 6) = -10*tau^3, (3, 7) = 15*tau^4, (3, 8) = -21*tau^5, (3, 9) = 28*tau^6, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1, (4, 5) = -4*tau, (4, 6) = 10*tau^2, (4, 7) = -20*tau^3, (4, 8) = 35*tau^4, (4, 9) = -56*tau^5, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1, (5, 6) = -5*tau, (5, 7) = 15*tau^2, (5, 8) = -35*tau^3, (5, 9) = 70*tau^4, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 1, (6, 7) = -6*tau, (6, 8) = 21*tau^2, (6, 9) = -56*tau^3, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 1, (7, 8) = -7*tau, (7, 9) = 28*tau^2, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 1, (8, 9) = -8*tau, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 1})

(2)

 


 

Download doMat.mw

is shown in the attached.

PS Have I seen this (or something very similar) before?

  restart;
  interface(rtablesize=10):
  interface(version);
  Physics:-Version();
  

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1)

  N:= 8:
  A:= Matrix
      ( N+1,
        N+1,
        (i,j)-> `if`
                 ( j>=i,
                   binomial(j-1, i-1)*(-tau)^(j-i),
                   0
                 )
      );

Matrix(9, 9, {(1, 1) = 1, (1, 2) = -tau, (1, 3) = tau^2, (1, 4) = -tau^3, (1, 5) = tau^4, (1, 6) = -tau^5, (1, 7) = tau^6, (1, 8) = -tau^7, (1, 9) = tau^8, (2, 1) = 0, (2, 2) = 1, (2, 3) = -2*tau, (2, 4) = 3*tau^2, (2, 5) = -4*tau^3, (2, 6) = 5*tau^4, (2, 7) = -6*tau^5, (2, 8) = 7*tau^6, (2, 9) = -8*tau^7, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = -3*tau, (3, 5) = 6*tau^2, (3, 6) = -10*tau^3, (3, 7) = 15*tau^4, (3, 8) = -21*tau^5, (3, 9) = 28*tau^6, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1, (4, 5) = -4*tau, (4, 6) = 10*tau^2, (4, 7) = -20*tau^3, (4, 8) = 35*tau^4, (4, 9) = -56*tau^5, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 1, (5, 6) = -5*tau, (5, 7) = 15*tau^2, (5, 8) = -35*tau^3, (5, 9) = 70*tau^4, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 1, (6, 7) = -6*tau, (6, 8) = 21*tau^2, (6, 9) = -56*tau^3, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 1, (7, 8) = -7*tau, (7, 9) = 28*tau^2, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 1, (8, 9) = -8*tau, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 1})

(2)

 


 

Download doMat.mw

as vv has already pointed out, with some functions you can use evalhf() and with other functions you can't.

Were you hoping for some kind of speed-up on the attached?


 

restart; interface(rtablesize = 10); hh := proc (cc::Array) local i; return add(pochhammer(cc[i], 2), i = 1 .. 1000) end proc; CodeTools:-Usage(evalf(hh(Array(1 .. 1000, [seq(h, h = 1 .. 1000)]))))

memory used=41.34MiB, alloc change=8.00MiB, cpu time=281.00ms, real time=286.00ms, gc time=0ns

 

334334000.0

(1)

 

NULL


 

Download evPoch.mw

 

The attached shows how to achieve "animations" by redefining the  display() function. Animating the "progress" of the geodesic is pretty simple. One thing I haven't been able to figure out is your logical basis for re-orienting the complete "tetrahedron+geodesic" plot, so I just made up a more-or-less random function for doing this latter aspect

The resulting animations display much better within Maple than they do on this site: it seems as if a couple mof the plot options aren't being correctly interpreted here - namely "transparency" so the geodesic appears and disappears, and "orientation" so the whole plot doesn't rotate.

You will have to take my word for the fact that both of these aspects will work within Maple if you download the worksheet (and the plots will just "look better" generally)

Geodesics on a tetrahedron

This worksheet provides tools for computing and

plotting geodesics on a tetrahedron.  The tetrahedron

need not be regular, although we take a regular tetrahedron

for the illustrations here.

restart;

with(plots):
interface(rtablesize=10):
interface(version);
Physics:-Version();

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1)

Vertices of a regular tetrahedron

pi[1], pi[2], pi[3], pi[4] := <0,0,0>, <1,0,0>, <1/2,sqrt(3)/2,0>, <1/2, sqrt(3)/6, sqrt(6)/3>:

Plot the tetrahedron with semi-transparent faces.

tetra := plottools:-tetrahedron(convert~([pi[1],pi[2],pi[3],pi[4]], list)):
T := display([
  display(tetra, style=line, color=black),
  display(tetra,color="PeachPuff", transparency=0.2, style=surface)
], scaling=constrained):

A face of the tetrahedron is specified through the indices u, v, w of its vertices.

E.g., u, v, w = 1, 4, 3corresponds to the face whose vertices are `&pi;__1`, `&pi;__4`, `&pi;__3`.

Along the edge connecting the vertices u and w we take a point P given by
P = a*`&pi;__u`+(1-a)*`&pi;__w`.

Along the edge connecting the vertices v and w we take a point Q given by

Q = b*`&pi;__v`+(1-b)*`&pi;__w`.

If we move from P to Q along the directed line segment PQ, and then

continue on the surface of the tetrahedron along a geodesic path, we

leave the previous face and enter a new face with vertices u__new, v__new, `w__new `

within which the geodesic traces a line segment QR.  The following proc

calculates and returns the indices u__new, v__new, w__new of the new vertices,

and the values of a and bcorresponding to the points Q and R on the

new face.

Later on, we construct a geodesic on the tetrahedron's surface through

repeated applications of this proc.

doit := proc(u, v, w, a, b)
  local u_new, v_new, w_new, a_new, b_new,
        fourth_vertex := ({1,2,3,4} minus {u,v,w})[],
        U := pi[u],
        V := pi[v],
        W := pi[w],
        P := a*U + (1-a)*W,
        Q := b*V + (1-b)*W,
        T := V+W-U,
        s1 := b*a/evalf((a-b)),
        t1 := -(b-1)*a/evalf(b);
  if t1 <= 1 then
    w_new := v;
    u_new := w;
    v_new := fourth_vertex;
    a_new := 1 - b;
    b_new := t1;
  else
    w_new := w;
    u_new := v;
    v_new := fourth_vertex;
    a_new := b;
    b_new := s1;
  end if;
  return u_new, v_new, w_new, a_new, b_new, Q;
end proc:

#
# Procedure to do required plot display
#
  doDisp:=proc(tet::function, P::list)
               uses plots:
               local j:
               return [ seq
                        ( display
                          ( [ tet,
                              pointplot3d( P[j],
                                           symbol=solidsphere,
                                           symbolsize=50,
                                           color=red
                                         ),
                              pointplot3d( P[1..j],
                                           style=line,
                                           color=blue
                                         )
                            ],
                            axes=none,
                            orientation=[160, 75, (360/numelems(P))*(j-1)]
                          ),
                          j=1..numelems(P)
                        )
                      ]:
        end proc:

Examples

 

In all of the following examples we begin with

u, v, w = 1, 2, 3, which means that the first line

segment of the geodesic lies in the `&pi;__1`, `&pi;__2`, `&pi;__3` plane

(which is the tetrahedron's base in our case).

 

Consider the vectors A = `&pi;__2`-`&pi;__1` and B = `&pi;__3`-`&pi;__1` as a basis

of the vectors in the plane of the base of the tetrahedron..

Let C = A*m+B*n.  For the first segment of the geodesic

we take a line within the base which parallels the vector C.

If the ratio m/n is rational, then the geodesic will be a closed

path, otherwise the geodesic curve will be dense in the

tetrahedron's surface.

 

The parameter a that appears in the computations below

may be assigned a value between 0 and 1. It specifies a

parallel displacement of the initial line segment

(which is parallel to the vector C as noted above.)  If a = 1,

the segment goes through the vertex `&pi;__1`.  If "a=0," it goes
through the vertex `&pi;__3`.  Don't change the parameter "b."

 

Example 1

 

We construct and plot 100 segments of the geodesic

but the result is a quadrilateral.  The plot is folded

over itself 25 times.

u:=1: v:=2: w:=3: m:=1: n:=0: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
#
# Replaced display command
#
  display( doDisp(T, [pts]), insequence=true);

 

Example 2

 

u:=1: v:=2: w:=3: m:=1: n:=1: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
#
# Replaced display command
#
  display( doDisp(T, [pts]), insequence=true);

 

Example 3

 

u:=1: v:=2: w:=3: m:=2: n:=3: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
#
# Replaced display command
#
  display( doDisp(T, [pts]), insequence=true);

 

Example 4

 

The ratio m/n is irrational, therefore the geodesic is dense

on the surface of the torus.

u:=1: v:=2: w:=3: m:=sqrt(2.0): n:=1: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
#
# Replaced display command
#
  display( doDisp(T, [pts]), insequence=true);

 

 


 

Download animGeodesic.mw

  1. I can't get really interested in "cosmetic" features
  2. Maple's typesetting is quite version-dependent, so what is shown in the attached may not reproduce in your (unspecified) version
  3. Specifying 'dot' and 'prime' notation (as far as I know) only works for univariate functions. For multivariate functions, 'compact' display is only available with subscripts.

Anyhow, see the attached for possibilities

  restart;

  interface(rtablesize=10):
  interface(version);
  Physics:-Version();

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 353, version installed in this computer: 350.`

(1)

#
# Maple on defaults
#
  diff(u(x),x)=2*diff(v(t),t);
  PDETools:-declare(u(x), y(t)):
  diff(u(x),x)=2*diff(v(t),t);

diff(u(x), x) = 2*(diff(v(t), t))

 

` u`(x)*`will now be displayed as`*u

 

` y`(t)*`will now be displayed as`*y

 

diff(u(x), x) = 2*(diff(v(t), t))

(2)

  restart;

#
# Juggle settings in Typesetting() package
#
# NB most of these only work for univariate
# functions - ie u(x) or u(t) but not for
# u(x,t)
#
  Typesetting:-Settings(dot=t):
  Typesetting:-Settings(usedot=true):
  Typesetting:-Settings(typesetdot=true):
  Typesetting:-Settings(prime=x):
  Typesetting:-Settings(typesetprime=true):
#
# Maple with typesetting adjustments
#
  diff(u(x),x)=2*diff(v(t),t);

diff(u(x), x) = 2*(diff(v(t), t))

(3)

 PDETools:-declare(u(x), y(t)):
 diff(u(x),x)=2*diff(v(t),t);

` u`(x)*`will now be displayed as`*u

 

` y`(t)*`will now be displayed as`*y

 

diff(u(x), x) = 2*(diff(v(t), t))

(4)

 

Download prettyDiff.mw

 

you can check the tensor as either

Weyl[alpha, beta, mu, nu];

or

Weyl[nonzero]

but not

Weyl[alpha, beta, mu, nu, nonzero]

See the attached

restart

with(Physics); Physics:-Version(); interface(rtablesize = 10); interface(version)

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 353, version installed in this computer: 350.`

 

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

(1)

Setup(mathematicalnotation = true)

[mathematicalnotation = true]

(2)

Setup(dimension = 3, coordinates = (X = [x1, x2, x3]), metric = 2*F6(X)*dx2*dx3+2*F5(X)*dx1*dx3+2*F4(X)*dx1*dx2+F1(X)*dx1^2+F2(X)*dx2^2+F3(X)*dx3^2)

`The dimension and signature of the tensor space are set to `[3, `- - +`]

 

`Default differentiation variables for d_, D_ and dAlembertian are:`*{X = (x1, x2, x3)}

 

`Systems of spacetime coordinates are:`*{X = (x1, x2, x3)}

 

_______________________________________________________

 

[coordinatesystems = {X}, dimension = 3, metric = {(1, 1) = F1(X), (1, 2) = F4(X), (1, 3) = F5(X), (2, 2) = F2(X), (2, 3) = F6(X), (3, 3) = F3(X)}]

(3)

g_[]

Physics:-g_[mu, nu] = Matrix(%id = 18446744074379560830)

(4)

Weyl[alpha, beta, mu, nu]; Weyl[nonzero]

0

 

Physics:-Weyl[alpha, beta, mu, nu] = {}

(5)

NULL

``

``

``

Download Weyl.mw

using the mehod shown in the attached.

The confusing(?) part of the definition in the referenced paper refers to 'N' being odd/even, but for the indexes in the matrices shown to make sense, then each matrix actually has to be N+1 X N+1

  restart;
  getMat:= N-> Matrix
               ( N+1,
                 N+1,
                 (i,j)-> `if`( type(j, odd),
                               `if`( `and`( type(i, odd),
                                            i>=j
                                          ),
                                      phi[ i-1, (i-j)/2 ],
                                      0
                                   ),
                               `if`( `and`( type(i, even),
                                            i>=j
                                          ),
                                      phi[ i-1, (i-j)/2 ],
                                      0
                                   )
                             )
               ):
  Z:= getMat(8); # N even
  Z:= getMat(9); # N odd

Matrix(9, 9, {(1, 1) = phi[0, 0], (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = phi[1, 0], (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = phi[2, 1], (3, 2) = 0, (3, 3) = phi[2, 0], (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (4, 1) = 0, (4, 2) = phi[3, 1], (4, 3) = 0, (4, 4) = phi[3, 0], (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (5, 1) = phi[4, 2], (5, 2) = 0, (5, 3) = phi[4, 1], (5, 4) = 0, (5, 5) = phi[4, 0], (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (5, 9) = 0, (6, 1) = 0, (6, 2) = phi[5, 2], (6, 3) = 0, (6, 4) = phi[5, 1], (6, 5) = 0, (6, 6) = phi[5, 0], (6, 7) = 0, (6, 8) = 0, (6, 9) = 0, (7, 1) = phi[6, 3], (7, 2) = 0, (7, 3) = phi[6, 2], (7, 4) = 0, (7, 5) = phi[6, 1], (7, 6) = 0, (7, 7) = phi[6, 0], (7, 8) = 0, (7, 9) = 0, (8, 1) = 0, (8, 2) = phi[7, 3], (8, 3) = 0, (8, 4) = phi[7, 2], (8, 5) = 0, (8, 6) = phi[7, 1], (8, 7) = 0, (8, 8) = phi[7, 0], (8, 9) = 0, (9, 1) = phi[8, 4], (9, 2) = 0, (9, 3) = phi[8, 3], (9, 4) = 0, (9, 5) = phi[8, 2], (9, 6) = 0, (9, 7) = phi[8, 1], (9, 8) = 0, (9, 9) = phi[8, 0]})

 

Matrix(%id = 18446744074378084710)

(1)

 

 


 

Download getMat.mw

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