tomleslie

5323 Reputation

15 Badges

9 years, 270 days

MaplePrimes Activity


These are answers submitted by tomleslie

Attached are two files.

The first file  (makeRepos.mw) creates the Maple repository and store it in the specified place. I have used full path specification for the location of the repository, so you may have to change the value of the variable 'reposit' to something more appropriate for your machine

I attached the relevant repository (ie foo.mla) to an email, then sent it to myself. I then saved the attachment to my desktop.

The second file (recLib.mw) shows how the recipient would modify his/her library path to execute the contents of 'foo.mla' (ie the procedure canSpell). Again I have used full path specification for the location of the saved email attachment, so the value of the variable 'myDesk' will have to be changed correspond to wherever the attachment was saved

#
# Create a temporary directory where I can store the
# Maple repository. I'm using the 'temp' directory
# on a windows machine.
#
# I'm using *full* paths to the directory, just to
# avoid any possible confusion/ambiguity.
#
# OP may want to change the definition of 'reposit'
# to something suitable for his/her machine
#
# If the directory already exists, remove it to ensure
# a nice "clean" start
#
  reposit := "C:/temp/reposit":
  if   FileTools:-Exists(reposit)
  then FileTools:-RemoveDirectory(reposit, recurse=true)
  fi;
#
# Now create the empty directory
#
  mkdir( reposit ):

#
# Create an empty repository called 'foo.mla' in the
# above directory.
#
# Again I'm using full paths, just to avoid any possible
# confusion/ambiguity
#
  march( 'create', cat( reposit, "/foo.mla"), 100 ):

#
# Check the contents of the repository "foo.mla". At this
# stage it should be empty!
#
  march('list', cat( reposit, "/foo.mla"));

[]

(1)

#
# Define the OP's procedure and save it to the repository
# created above
#
  canSpell := proc(w)
                   option encrypted;
                   local blocks, i, j, word, found, N;
                   blocks := Array([{"B", "O"}, {"X", "K"}, {"D", "Q"}, {"C", "P"}, {"N", "A"},
                                    {"J", "W"}, {"H", "U"}, {"V", "I"}, {"A", "N"}, {"O", "B"},
                                    {"G", "T"}, {"R", "E"}, {"T", "G"}, {"Q", "D"}, {"F", "S"},
                                    {"E", "R"}, {"F", "S"}, {"L", "Y"}, {"P", "C"}, {"Z", "M"}]);
                   word := StringTools:-UpperCase(convert(w, string));
                   N := numelems(blocks);
                   for i to length(word) do
                       found := false;
                       for j to N do
                           if word[i] in blocks[j] then
                              blocks[j] := blocks[N];
                              N := N-1;
                              found := true;
                              break;
                           end if;
                       end do;
                       if not found then
                          return false;
                       end if;
                   end do;
                   return true;
             end proc:

  savelib(canSpell, cat( reposit, "/foo.mla")):

#
# Recheck the contents of the repository
#
# Everything *seems* to be working
#
  march('list', cat( reposit, "/foo.mla"));
  

[["canSpell.m", [2019, 9, 19, 11, 32, 46], 41984, 1814]]

(2)

 

Download makeRepos.mw


 

  restart;

#
# After creating the repository, I attached it to an email
# and sent it to myself. I then saved the attachment to my
# "Desktop", defined below.
#
# The value of 'myDesk' below will need to be changed by the
# recipient of OP's email to correspond to where (s)he saved
# the email attachment
#
  myDesk:="C:/Users/TomLeslie/Desktop":

#
# Check my default library path
#
  libname;
#
# Temporarily redefine the libname to include the "foo.mla"
# library on my desktop
#
  libname:=libname, cat( myDesk, "/foo.mla");

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib", "C:\Program Files\Maple 2019\lib", "C:\Users\TomLeslie\maple\toolbox\OrthogonalExpansions\lib", "C:\Users\TomLeslie\maple\toolbox\Syrup\lib"

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib", "C:\Program Files\Maple 2019\lib", "C:\Users\TomLeslie\maple\toolbox\OrthogonalExpansions\lib", "C:\Users\TomLeslie\maple\toolbox\Syrup\lib", "C:/Users/TomLeslie/Desktop/foo.mla"

(1)

#
# Check the contents of the added repository
#
  march('list', cat( myDesk, "/foo.mla") );

[["canSpell.m", [2019, 9, 19, 10, 41, 59], 41984, 1814]]

(2)

#
# Check encryption
#
  eval(canSpell);

proc (w) local blocks, i, j, word, found, N; option encrypted; "2f0b0404264baa44a4bc758aa04e5b18e763f64743733b1443d884fbae3db43aa9f7602fc86ade456f648e4965c66765bc9ae483a741f33e7906cb36c2c9f8d8571c607892e334a524080dc68f96c7ee03127724368a0b44518f210cc61a58ef47668d7415828c6e4de386370634e04754d32a723d154802511873d5b144d08ba5c8d544356a5126d8d39e047de4fb7251107155148450fca5c8d5443d6a51265be4a391eb29d40e51197165d872d08ba5c8d5443a6a512690f1bf332a8a69ef511a71556c8450fc402c48ce1744215c81a7ede7e8baefb3a5c8d544396a5126208cee0573051de2ff10f544e08e25c34119983289241252bdbc456f5124c54684d0a607a99430dae81d6220cc795ad05e54c1fe2ab47fa3d4b4515eb853703920a8ac666e690a744934462441290b6e62016dc6f0f1d3dc45b49f2d356d666ebcf3c846116cc4ef57d407aeef5ccc29c3c4a776c0780b78a86eabbfffc16dc83ceb6e4776963171fe468b8fde1fcfbbe01d9f706c057c45b2803d6931787b36c93468625240bced45e82c676e653427f554f621d038706303149886f5921caac96f01288b3e65433f340647989334611766b05d3938b5f3e7454c42c6377bade609147920fb7524357aa946133c235b6c5f0da1705c06234890d7203931aee0dcbdb9463b07c955262c3767463e376773c1580fd747af44f886db6e4e5d4e213c1858365a9155f3464d8865464f286496b060ab4849d13be286ee895dccbb34b5702271feedf9252c51763413f9a02cb2a23115a91b68458c65f69193db79837264567113463261dccaecb83b95a6523e225ba3287dc64429c643319c1f94358356eb6f54d4492fd82582a610073f20c13b2a5f28a877ae735d5671f3463261ba0717e67c56741447c81c48824725a3f0a2f339af8153f5d669f00d697dcfad76c0cd0837e068ebd6896ad76e2ee288736556710b4632616a39f97a47d8f655c1225b1429103e4481e369e6d4a3aecab2458684dfc0612dd82582a610073f20a9100e6b231e74197264567113463261af062879482fde19" end proc

(3)

#
# Run the OP's procedure
#
  canSpell("good");

true

(4)

 


 

Download recLib.mw

Do all calculations with complex numbers

Just need a a couple of functions which convert inputs from "phasor" to complex and then converts output from complex back to phasor, as in the following

#
# Define a couple of simple functions which
#
# 1) Convert from "phasor" representation to
#    complex numbers
# 2) Convert from complex representation to
#    "phasor" representation
#
  fromP:= (r, theta)-> r*exp(I*Pi*theta/180);
  toP:= x-> evalf([abs(x), (180/Pi)*argument(x)]);
#
# Define "phasors" (and convert to complex
# reprentation)
#
  U__LI:= fromP(230, 0):
  I__LI:= fromP(20, -30):
  Z__L:=  fromP(0.097, 7.2):
  I__N:=  fromP(40,-120):
#
# Perform calculation, and convert to "phasor"
# representation
#
   ans:= toP(U__LI - I__LI * Z__L + I__N * Z__L);

proc (r, theta) options operator, arrow; r*exp(((1/180)*I)*Pi*theta) end proc

 

proc (x) options operator, arrow; evalf([abs(x), 180*argument(x)/Pi]) end proc

 

[226.7256261, -.7139358893]

(1)

 


 

Download phasors.mw

In the attached

  1. I removed some "comment" areas which appeared to be a mixture of "comment" and "executable code" - always tricky to tell the difference if you insist on using 2-D input
  2. I didn't like the syntax for the boundary condition D(u(0, t)) = 0, Slightly surprised that Maple didn't object to this, but I changed it to something which is definitely syntactically correct.

Even then I still get an error that

Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0): matrix is singular

Now I can't be certain why this error message is occurring, but since the problem seems to occur at (x,t)=(0,0), one thing that bothers me is the BCs/ICs at (0,0).

You have the initial condition

u(x, 0) = 0.1*sin(2*Pi*x)

which must be true for all 'x' - so the above requires that u(0,0)=0.

However you also have the boundary condtion

u(0, t) = 0.5*cos(2*Pi*t)

which is true for all 't', so this requires that u(0,0)=0.5 in conflict with initial condition above

A similar conflict occurs with the BCs/ICs for v(0,0): the relevant initial condition has v(0,0)=0, wheres the relevant boundary condition has v(0,0)=0.2

You might want to give some thought to these boundary/initial condition conflicts

m := 1; hBar := 1

pdeu := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -0.1e-1

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -0.1e-1

(1)

pdev := diff(v(x, t), t)+u(x, t)*(diff(v(x, t), x))-hBar*(diff(u(x, t), `$`(x, 2)))/(2*m)+v(x, t)*(diff(u(x, t), x))/m = -0.2e-1

diff(v(x, t), t)+u(x, t)*(diff(v(x, t), x))-(1/2)*(diff(diff(u(x, t), x), x))+v(x, t)*(diff(u(x, t), x)) = -0.2e-1

(2)

IC := {u(x, 0) = .1*sin(2*Pi*x), v(x, 0) = .2*sin((1/2)*Pi*x)}

{u(x, 0) = .1*sin(2*Pi*x), v(x, 0) = .2*sin((1/2)*Pi*x)}

(3)

BC := {u(0, t) = .5*cos(2*Pi*t), v(0, t) = .2*cos(2*Pi*t), (D[1](u))(0, t) = 0}

{u(0, t) = .5*cos(2*Pi*t), v(0, t) = .2*cos(2*Pi*t), (D[1](u))(0, t) = 0}

(4)

``

IBC := `union`(IC, BC)

{u(0, t) = .5*cos(2*Pi*t), u(x, 0) = .1*sin(2*Pi*x), v(0, t) = .2*cos(2*Pi*t), v(x, 0) = .2*sin((1/2)*Pi*x), (D[1](u))(0, t) = 0}

(5)

pds := pdsolve({pdeu, pdev}, IBC, numeric, time = t, range = 0 .. 1)

_m747256992

(6)

p1 := pds:-plot(t = 10, numpoints = 1000, color = blue)

Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
matrix is singular

 

NULL

Download pdeProb2.mw

In the atached,  removed a few things you probably(?) don't need.

Everything seems to be working correctly

What exactly do you want to do?

restart

with(plots)

L := 1; alpha := 1; epsilon := .1; `𝒳` := x^2

Inverse := diff(C(x, t), t) = -epsilon*(diff(C(x, t), x, x, x, x))-alpha*(diff(C(x, t), x, x))

IBCI := C(x, 0) = `𝒳`, C(0, t) = 0, C(L, t) = 0, (D[1, 1](C))(0, t) = 0, (D[1, 1](C))(L, t) = 0

pdsI := pdsolve(Inverse, [IBCI], numeric); pdsI:-plot3d(C(x, t), t = 0 .. 10, x = 0 .. 1); pdsI:-plot(t = 10, x = 0 .. 1); pdsI:-plot(x = .25, t = 0 .. 10)

_m696589984

 

 

 

 

``


 

Download pdeProb.mw

when performing multiple substitutions, like x=..., t=.., then these have to grouped in a list, as in [x=.., t=..] as in the attached

For future reference, when posting here, please use the big green up-arrow in the Mapleprimes toolbar to upload code. 'Pictures' of code are not editable, not executable, and so not very useful

  restart:

  f:=7*exp(-3*t)-(3*x):
  s:=0.05:
  T:=0.1:
  n:=round(T/s):
  t0:=0:
  x[0]:=6:
  for i from 0 by 1 to n-1 do
    #
    # Note the square brackets enclosing
    # the multiple substitutions
    #
    # the simplify() 'wrapper' is just so that
    # exp(0) gets evaluated to a floating point
    # which isn't really necessary at this stage
    #
      x[i+1]:=simplify(x[i]+s*subs( [x=x[0], t=t0+i*s], f)):
  od;

5.45

 

4.851247792

(1)

restart

f := 7*exp(-3*t)-3*x
s := 0.5e-1
T := .1
n := round(T/s)
t0 := 0

x[0] := 6

for i from 0 to n-1 do x[i+1] := simplify(x[i]+s*subs([x = x[0], t = i*s+t0], f)) end do

5.45

 

4.851247792

(2)

 

Download doSub.mw

Try adding more-or-less plausibe constraints to the values of the solution variables - as in the attached

NB I don't know enough about fluids to determine whether the constraints in the attached are plausible or not - I just guessed. You may be able to "relax" these constraints quite a bit and still get a solution. As a general rule, any constraint you can add to an fsolve() problem (even something as simple as requiring all variables to be positive) increase the chances of obtaining a solution

Anyhow, for what it is worth, the second fsolve() command in the attached "works" in Maple 2019

  restart;
  interface(version);

`Standard Worksheet Interface, Maple 2019.1, Windows 7, May 21 2019 Build ID 1399874`

(1)

  with(ThermophysicalData);
  with(ThermophysicalData[CoolProp]);
  fluid := "R717"; T__in := 310; p__in := 11*10^5; v__in := 10;  A := 0.005;
  h__in := Property(massspecificenthalpy, T = T__in, P = p__in, fluid);
  rho__in := Property(density, T = T__in, P = p__in, fluid);
  m := A*v__in*rho__in;
  p__out := 2*10^5;
  eq1 := h__out = Property("massspecificenthalpy", "temperature" = T__out, "P" = p__out, fluid);
  eq2 := rho__out = Property("D", "temperature" = T__out, "P" = p__out, fluid);
  eq3 := h__in + v__in^2/2 = h__out + v__out^2/2;
  eq4 := m = A*v__out*rho__out;
  res := fsolve({eq1, eq2, eq3, eq4});
#
# Try adding plausible(?) constraints on the ranges
# of the variables
#
  res := fsolve( {eq1, eq2, eq3, eq4},
                 {T__out=200..500, h__out=0..5e06, v__out=0..100, rho__out=0..5}
               );

[Chemicals, CoolProp, PHTChart, Property, PsychrometricChart, TemperatureEntropyChart]

 

[HAPropsSI, PhaseSI, Property, Props1SI, PropsSI]

 

"R717"

 

310

 

1100000

 

10

 

0.5e-2

 

1655590.94730295590

 

8.15132757401520713

 

.4075663787

 

200000

 

h__out = ThermophysicalData:-CoolProp:-Property("massspecificenthalpy", "temperature" = T__out, "P" = 200000, "R717")

 

rho__out = ThermophysicalData:-CoolProp:-Property("D", "temperature" = T__out, "P" = 200000, "R717")

 

1655640.947 = h__out+(1/2)*v__out^2

 

.4075663787 = 0.5e-2*v__out*rho__out

 

fsolve({.4075663787 = 0.5e-2*v__out*rho__out, 1655640.947 = h__out+(1/2)*v__out^2, h__out = ThermophysicalData:-CoolProp:-Property("massspecificenthalpy", "temperature" = T__out, "P" = 200000, "R717"), rho__out = ThermophysicalData:-CoolProp:-Property("D", "temperature" = T__out, "P" = 200000, "R717")}, {T__out, h__out, v__out, rho__out})

 

{T__out = 285.0179004, h__out = 1654113.289, v__out = 55.27490598, rho__out = 1.474688637}

(2)

 

Download tProp.mw

 

 

something like the attached - although I don't think you gain much (if anything) when compared to using 'D-operator' notation


 

  restart;

#
# Define function
#
  doDiff:= (p,q,r)-> diff(p(q), q$r);
#
# A couple of examples
#
  doDiff( f, x, 3);
  doDiff( g, y, 2);
  doDiff( h, z, n);

proc (p, q, r) options operator, arrow; diff(p(q), `$`(q, r)) end proc

 

diff(diff(diff(f(x), x), x), x)

 

diff(diff(g(y), y), y)

 

diff(h(z), [`$`(z, n)])

(1)

#
# But is this any better that using
# D-operator notation?
#
  D[1$3](f)(x);
  D[1$2](g)(y);
  D[1$n](h)(z);
  convert(%%%, diff);
  convert(%%%, diff);
  convert(%%%, diff);

((D@@3)(f))(x)

 

((D@@2)(g))(y)

 

((D@@n)(h))(z)

 

diff(diff(diff(f(x), x), x), x)

 

diff(diff(g(y), y), y)

 

diff(h(z), [`$`(z, n)])

(2)

 

 


 

Download doDiff.mw

to use the command

MmaTranslator[FromMmaNotebook]()

which according to its help page

The FromMmaNotebook(Mma_notebook_filename) and the convert(Mma_notebook_filename, FromMmaNotebook) commands translate a Mathematica notebook to a Maple worksheet and saves the results to disk. These commands enable Mathematica users to automatically translate their Mathematica notebooks to Maple worksheets.

 

but can't find it?!

So far as I can tell, an analytic solution is not possible - so one is forced to look for a numeric solution

In a numeric solution, it is not possible to use "infinity" as one of the boundaries. The "conventional" approach is to use an appropriately "large" number as a "proxy" for infinity.

Even doing this, I had various problems producing a numeric solution which "worked" over a range of values for 'delta', which resulted in the selection of a few specific options for the dsolve(....., numeric) command.

The attached *appears* to be OK for the required range of delta-values, including delta=0

Anyhow for what it is worth, check the attached.

NB I couldn't check this in Maple 17, because the earliest version I have loaded is Maple 18

  restart;
  with(plots):
  delta:=1:
  inf:=500*delta+1;
  sol:= dsolve( [ 2*diff(y(x),x,x,x)+y(x)*diff(y(x),x,x)=0,
                  y(0)=0,
                  D(y)(inf)=1,
                  D(y)(0)=delta*D[1,1](y)(0)
                ],
                numeric,
                approxsoln=[y(x)=x],
                maxmesh=4096,
                method=bvp[midrich]
              ):
  odeplot( sol, [x, y(x)], x=0..1);
  odeplot( sol, [x, y(x)], x=0..inf)

501

 

 

 

 


 

Download BVPprob.mw

For linear systems, one can use commands from the LinearAlgebra package to examine whether the system is consistent or underdetermined etc

As in the attached

  restart:
#
# Make sure enough digits are used for calculation
# but restrict how many are displayed, just for
# readibility
#
  Digits:=50:
  interface(displayprecision=4):
#
# Since system is 'linear' probably(?) better to use
# LinearAlgebra:-LinearSolve() rather than the more
# "straightforward" fsolve(). so load the relevant
# package
#
  with(LinearAlgebra):
#
# OP's system
#
  C[0]:= 3.19153824321146142351956847947*tau[1]-19.1492294592687685411174108768*tau[2]+111.703838512401149823184896781*tau[3]+3.19153824321146142351956847947*tau[4]-44.6815354049604599292739587124*tau[5]+622.349957426234977586315853494*tau[6]:
  C[1]:= 51.0646118913833827763130956714*tau[2]-612.775342696600593315757148056*tau[3]+51.0646118913833827763130956714*tau[5]-1429.80913295873471773676667880*tau[6]:
  C[2]:= -1.06073680388443795908856507616+3.19153824321146142351956847947*tau[1]+53.1609155734306093706448370717*tau[2]+1672.89412862088744108725223170*tau[3]+3.19153824321146142351956847947*tau[4]+27.6286096277389179824882892361*tau[5]+1026.57792701153122226218722129*tau[6]:
  C[3]:= -1.08847004231036963538035920033+3.19153824321146142351956847947*tau[1]+62.6399144226357196540662623767*tau[2]+2040.52109049201342887896297462*tau[3]+3.19153824321146142351956847947*tau[4]+37.1076084769440282659097145411*tau[5]+1242.54090729537544551915515930*tau[6]:
  C[4]:= -1.05523181556926815105314303389+3.19153824321146142351956847947*tau[1]+72.7671212023804312453829273862*tau[2]+2472.93216226733267613216245895*tau[3]+3.19153824321146142351956847947*tau[4]+47.2348152566887398572263795506*tau[5]+1512.91667059477930731128800348*tau[6]:
  C[5]:= -.922876006485286011069063957991+3.19153824321146142351956847947*tau[1]+82.9822841707707093164204255644*tau[2]+2971.36790137532483139495115633*tau[3]+3.19153824321146142351956847947*tau[4]+57.4499782250790179282638777288*tau[5]+1847.90980220852701343747673000*tau[6]:
#
# Convert to a matrix system and attempt a solution
#
  A,b:= GenerateMatrix
        ( [seq( C[j], j=0..5)],
          [seq( tau[j], j=1..6)]
        );
  LinearSolve(A,b);

A, b := Matrix(6, 6, {(1, 1) = 3.19153824321146142351956847947, (1, 2) = -19.1492294592687685411174108768, (1, 3) = 111.703838512401149823184896781, (1, 4) = 3.19153824321146142351956847947, (1, 5) = -44.6815354049604599292739587124, (1, 6) = 622.349957426234977586315853494, (2, 1) = 0, (2, 2) = 51.0646118913833827763130956714, (2, 3) = -612.775342696600593315757148056, (2, 4) = 0, (2, 5) = 51.0646118913833827763130956714, (2, 6) = -1429.80913295873471773676667880, (3, 1) = 3.19153824321146142351956847947, (3, 2) = 53.1609155734306093706448370717, (3, 3) = 1672.89412862088744108725223170, (3, 4) = 3.19153824321146142351956847947, (3, 5) = 27.6286096277389179824882892361, (3, 6) = 1026.57792701153122226218722129, (4, 1) = 3.19153824321146142351956847947, (4, 2) = 62.6399144226357196540662623767, (4, 3) = 2040.52109049201342887896297462, (4, 4) = 3.19153824321146142351956847947, (4, 5) = 37.1076084769440282659097145411, (4, 6) = 1242.54090729537544551915515930, (5, 1) = 3.19153824321146142351956847947, (5, 2) = 72.7671212023804312453829273862, (5, 3) = 2472.93216226733267613216245895, (5, 4) = 3.19153824321146142351956847947, (5, 5) = 47.2348152566887398572263795506, (5, 6) = 1512.91667059477930731128800348, (6, 1) = 3.19153824321146142351956847947, (6, 2) = 82.9822841707707093164204255644, (6, 3) = 2971.36790137532483139495115633, (6, 4) = 3.19153824321146142351956847947, (6, 5) = 57.4499782250790179282638777288, (6, 6) = 1847.90980220852701343747673000}), Vector(6, {(1) = 0, (2) = 0, (3) = 1.06073680388443795908856507616, (4) = 1.08847004231036963538035920033, (5) = 1.05523181556926815105314303389, (6) = .922876006485286011069063957991})

 

Error, (in LinearAlgebra:-BackwardSubstitute) inconsistent system

 

#
# Idle curiosity - check the reduced row echelon form
# of the augmented matrix. Three observations
#
# 1. the fact that row 6 is all zero implies that the
#    system is underdetermined - so any solution would
#    require a free parameter, however
# 2. row 5 is equivalent to the statement 0=1, so no
#    solution is possible
# 3. The fact that all of the coefficients in the RRF
#    *seem* to be *exact* makes one wonder exactly how
#    the original equation system was generated!
#
  M:=ReducedRowEchelonForm(<A|b>);
#
# Regenerate the system of equations from the reduced
# row echelon form. Check the fifth and sixth entry
# in the output list, which confirms the statements
# (1) and (2) above
#
  GenerateEquations(M, [seq( tau[j], j=1..6)]);

Matrix(6, 7, {(1, 1) = 1., (1, 2) = 0., (1, 3) = 0., (1, 4) = 1.0000000000000000000000000000000000000000000000000, (1, 5) = -7.9999999999999999999999999999498674345073799899517, (1, 6) = 0., (1, 7) = 0., (2, 1) = 0., (2, 2) = 1., (2, 3) = 0., (2, 4) = 0., (2, 5) = 1.0000000000000000000000000000000000000000000000000, (2, 6) = 0., (2, 7) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 1., (3, 4) = -0., (3, 5) = -0., (3, 6) = -0., (3, 7) = -0., (4, 1) = 0., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0., (4, 5) = 0., (4, 6) = 1., (4, 7) = 0., (5, 1) = 0., (5, 2) = 0., (5, 3) = 0., (5, 4) = 0., (5, 5) = 0., (5, 6) = 0., (5, 7) = 1., (6, 1) = 0., (6, 2) = 0., (6, 3) = 0., (6, 4) = 0., (6, 5) = 0., (6, 6) = 0., (6, 7) = 0.})

 

[1.*tau[1]+1.0000000000000000000000000000000000000000000000000*tau[4]-7.9999999999999999999999999999498674345073799899517*tau[5] = 0., 1.*tau[2]+1.0000000000000000000000000000000000000000000000000*tau[5] = 0., 1.*tau[3] = -0., 1.*tau[6] = 0., 0. = 1., 0. = 0.]

(1)

 

 

Download linSol.mw

 is shown in the attached (and I am sure there are many other ways!)

  restart;
#
# Implement the sieve of Eratosthenes. This will
# return all prime numbers up to the supplied
# argument (reasonably efficiently)
#
# EG on my machine all primes up to 1000000 can
# be generated in less than 1sec
#
  sieve := proc( Nval::integer)
                 local vec, k;
                 description "The Sieve of Eratosthenes";
                 uses ArrayTools, LinearAlgebra:
               #
               # Initialise a vector with all ones
               #
                 vec:= Vector[row]( Nval, 1 ):
               #
               # set the first entry in the vector to zero
               #        
                 vec[1]:= 0:
               #
               # Set all even entries with index >=4, to zero
               # (eliminates all the even indices which can't
               # be prime)
               #
                 Fill(0, vec, 3, 2);
               #
               # For each odd number, say k,  eliminate all
               # of its odd multiples (other than 1*k)
               #
                 seq( Fill
                      ( 0, vec, k^2-1, k ),
                      k=3..floor( sqrt(Nval) ),
                      2
                    ):
               #
               # return the indices for all non-zero values
               #
                 return Transpose
                        ( SearchArray
                          ( vec )
                        );
           end proc:
#
# Generate all the primes between 1 and N, Change
# the following assignment as necessary
#
  N:=1000000:
  prims:=sieve(N):
#
# Generate a random number between 1 and the
# number of primes in the vector 'prims'
#
  randomize():
  r:= rand( 1..op(1, prims)):

#
# Now generate (say) a hundred "random" prime
# numbers between 1 and N (defined above)
#
  seq( prims[ r() ], j=1..100);

575119, 538457, 296669, 748637, 180311, 486757, 223753, 52673, 854083, 257, 651769, 995173, 499253, 878641, 479287, 910603, 294923, 318863, 373151, 704251, 174487, 151471, 620183, 491527, 603391, 184241, 174859, 148763, 505339, 701653, 870031, 141269, 696433, 390107, 927833, 42061, 102533, 678779, 141679, 183377, 100537, 982687, 54091, 227363, 681061, 678581, 372013, 533909, 461861, 276883, 56591, 346039, 102233, 131743, 794413, 639371, 712493, 46327, 874763, 273719, 15307, 632743, 151483, 515639, 941209, 580913, 695677, 599303, 56311, 739217, 560137, 722411, 603613, 994067, 377623, 19937, 70241, 488321, 135449, 525949, 397939, 418909, 373, 125933, 909281, 465089, 840589, 201829, 268607, 104851, 559877, 967507, 44771, 829537, 24121, 911173, 215959, 625267, 672209, 958501

(1)

 

Download primRand.mw

Put anything you want/need in your Maple initialization file. Seee the help at

?Create Maple Initialization file

You will already have such a file - although there is probably(?) noting in it

The differential transformation method is essentially a power series solution (some authors appear to think that it is actually a Taylor series method in disguise?!). It would seem to be appropriate for an initial-value problem, where the expansion point of the power series is the initial vaue.

For a BVP with defined at two (boundary) points, about whihc of these points would the power series expansion take place? Hence the DTM would seem to be completely inappropriate for general BVPs.

For your example, a DTM method appears not only inappropriate, but unnecessary. One only has to decide how to handle the boundary condition at" infinity". The conventional approach is to choose the value of the second boundary to be a "large" number, (eg try 10, 100, 1000) as a substitute for "infinity" and see if the solution appears to settle down to something "reasonable".

The attached code appears to work for most values of the parameter delta and doesn't change form much if one sets the upper limit to 10,100,1000. The selected options in the dsolve() command were necessary to ensure a solution for the variious values of 'delta' and the selected upper limit

delta:=-5:
inf:=1000;
ode:=2*diff(y(x),x,x,x)+y(x)*diff(y(x),x,x)=0:
ics:= y(0)=0, D(y)(inf)=1, D(y)(0)=delta*D[1,1](y)(0):
sol:=dsolve( [ode, ics], numeric, approxsoln=[y(x)=x], maxmesh=4096,method=bvp[midrich]);
odeplot( sol, [x, y(x)], x=0..1);
odeplot( sol, [x, y(x)], x=0..inf)

1000