tomleslie

5189 Reputation

15 Badges

9 years, 244 days

MaplePrimes Activity


These are answers submitted by tomleslie

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