tomleslie

6237 Reputation

17 Badges

10 years, 106 days

MaplePrimes Activity


These are replies submitted by tomleslie

@Carl Love 

  1. "Floating point subscripts" I can make this behaviour "go away" if I export the OP's original as a .mpl file and then read() the .mpl file in a new worksheet. I only tried this on the off-chance that there might be some "non-printing" character in the original, which would either break the export process, or be identifiable in an external editor. However the export process worked as normal: and the .mpl file showed no sign of non-ascii characters in an external editor  - so pretty inconclusive.
  2. Why does the revised worksheet  oddEqs.mw  I posted obtain a non-zero result? This seems to depend on the precise  "preprocessing" of the equations prior to the solve() command. OP uses combine(expand(...)) where as my original worksheet omitted the combine(..) step. Since I was using solve(), to provide an exact solution, I thought the combine() step was redundant, although as the attached shows - it does make a significant difference to the solution(s) obtained. I then thought this might be  some floating-point-error phenomenon (might still be??), but changing the 'Digits' setting doesn't seem to make a difference - without the combine() step, a non-zero answer is obtained: with the combine() step only a zero answer is obtained. See the attached.

  restart:

  odeA1:=   0.0001995000000*a__1*omega^2*cos(omega*t)
          + 0.0001995000000*a__3*omega^2*cos(3*omega*t)
          + 0.004987500000*a__5*omega^2*cos(5*omega*t)
          - 93.86*a__1*cos(omega*t)
          - 10.42888889*a__3*cos(3*omega*t)
          - 93.86*a__5*cos(5*omega*t)
          + 1.633844046*10^18*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))^2*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9 + a__5*cos(5*omega*t))
          + 0.02250237053*a__1*omega*sin(omega*t)
          + 0.007500790177*a__3*omega*sin(3*omega*t)
          + 0.1125118526*a__5*omega*sin(5*omega*t):
  odeA2:=   A__1*omega*sin(omega*t)
          + A__3*omega*sin(3*omega*t)/3
          + 5*A__5*omega*sin(5*omega*t)
          - 1.499250375*10^(-7)*(A__1*cos(omega*t)
          + A__3*cos(3*omega*t)/9
          + A__5*cos(5*omega*t))*(2.110712380*10^10
          - 1.633844050*10^18*(a__1*cos(omega*t)
          + a__3*cos(3*omega*t)/9
          + a__5*cos(5*omega*t))^2)
          + 1.499250375*10^(-7):
#
# With the 'combine()' step
#
  sol1:= eval
         ( solve
           ( [ seq
               ( coeff~
                 ( combine~
                   ( expand~
                     ( [odeA1,odeA2]
                     )
                   ),
                   cos(j*omega*t)
                 )[],
                 j=1..15,2
               )
             ],
             [ a__1, a__3, a__5, A__1, A__3, A__5]
           )[],
          omega=485.0135000
        );
#
# Without the 'combine()' step
#
  sol2:= eval
         ( solve
           ( [ seq
               ( coeff~
                 ( expand~
                   ( [odeA1,odeA2]
                   ),
                   cos(j*omega*t)
                 )[],
                 j=1..15,2
               )
             ],
             [ a__1, a__3, a__5, A__1, A__3, A__5]
           )[],
           omega=485.0135000
         );

[a__1 = 0., a__3 = 0., a__5 = 0., A__1 = 0., A__3 = 0., A__5 = 0.]

 

[a__1 = -2.333333333*a__3+115.0000000*a__5, a__3 = a__3, a__5 = a__5, A__1 = A__1, A__3 = 2.999999999*A__1+15.*A__5, A__5 = A__5]

(1)

 

NULL


Download oddEqs2.mw
 

 

@erik10 

my original answer.

  1. No-one apart from you has reported a problem opening .mw files in Maple 2020 in Windows 10
  2. Quite a lot of people seem to have a problem with broken file associations in Windows 10. This is a Windows 10 problem
  3. So somehow your Windows installation won't allow/accept/retain/setup a file association between Maple 2020 and .mw files
  4. Can only suggest you Google "Problems with Windows 10 file associations" and start ploughing through the many suggestions

@raj2018 

and by "complex" I don't mean "complicated", but (literally) it generates complex numbers. Now Maple can solve IVP problems whihc generate complex solutions - but pretty much everything gets more complicated, very quickly.

A trivial instance of the complication - since solutions of A(x) are complex, what does one choose as a suitable initial condition, since A(0) is now allowed (required?) to be a complex number?

The attached shows a couple of solution cases with A(x) complex, for real initial conditions.

Before I go any further with this, a number of questions need answering, such as

  1. Are you expecting A(x) to be complex? If not, why not?
  2. Since A(x) appears to be generally complex, what is an acceptable value for the (possibly complex) initial value A(0). Whether (or not) numerical solutions can be obtained, seems to be very dependent on this initial condition
  3. If an ODE is providing solutions that you really (really) do not expect, this *usualll* means
    1. The problem is incorrectly set up - maybe some typo somewhere - just crawl all over everything
    2. You are (without realising) making some kind of assumptions - which of course Maple will not do. So think long and hard about why you are assuming certain kinds of solutions

  restart:
  q1:=0.4: q:=0.4: s1:=0: mu:=2: mu1:=0.7:
  mu3:=0.2: s:=0.5: M:=1.9109: mu2:=1+mu3-mu1:

  V:=unapply( M^(2)-M*(M^(2)+2*A)^(1/(2))+mu1*(M^(2)*mu+s1)-mu1*((M^(2)*mu)/(2))^(1/(2))*(((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(1/(2))+(4*M^(2)*mu*s1)*((M^(2)*mu+3*s1-2*A)+sqrt((M^(2)*mu+3*s1-2*A)^(2)-12*M^(2)*mu*s1))^(-3/(2)))-(2*mu2)/((3*q-1))*((1-(q-1)*A)^((3*q-1)/(2*q-2))-1)-(2*mu3)/(s*(3*q1-1))*((1+s*(q1-1)*A)^((3*q1-1)/(2*q1-2))-1),A):

#  ode:=diff(A(x),x)^2+V(A(x))=0;
#
# Manually rewrite ODE selecting one sign of the
# square root
#
  ode:=diff(A(x),x)=sqrt(V(A(x))):
  odeSol:=dsolve([ode, A(0)=1], numeric, complex=true):
#
# Plot the real and imaginary parts
#
  plots:-odeplot( odeSol, [ [x, Re(A(x))],[x, Im(A(x))] ], x=0..10, color=[red,blue]);

 

 

#
# Now use the other sign of the square roo
#
  ode:=diff(A(x),x)=-sqrt(V(A(x))):
  odeSol:=dsolve([ode, A(0)=-1], numeric, complex=true):
#
# Plot the real and imaginary parts
#
  plots:-odeplot( odeSol, [ [x, Re(A(x))],[x, Im(A(x))] ], x=-10..0, color=[red,blue]);

 

 

 

Download odeComplex.mw

@AHSAN 

I have written three (very simple) rootfinding procedures, using three conceptually different iteration processes, namely

  1. regula falsi
  2. bisection
  3. newton's method

  restart;
#
# Procedure which computes 'root' using regula falsi,
# returns root, residual, and number of iterations
#
  regFals:= proc( g, X__A, X__B, err)
                  local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
                  uses Student[Precalculus]:
                #
                # Iterate until tolerance is achieved
                #
                  for j from 1 do
                      iters:=iters+1:
                    #
                    # Get the x-intercept of the line between
                    # the two points  given by {X__1, f(X__1)]
                    # and [X__2, f(X__2)]
                    #
                      p:=evalf( [ Line
                                  ( [X__1, g(X__1)],
                                    [X__2, g(X__2)]
                                  )
                                ][4]
                              );
                    #
                    # Exit if required tolerance has be achieved
                    #
                      if   abs(g(p))<err
                      then break
                    #
                    # Otherwise update the range end-points
                    #
                      elif g(p)>0
                      then X__2:=p:
                      else X__1:=p;
                      fi;
                 od:
                 return root=p, residual=g(p), numberOfIterations=iters;
            end proc:
#
# Procedure which computes root using the bisection method
#
  bisection:= proc( g, X__A, X__B, err )
                    local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
                  #
                  # Iterate until tolerance is achieved
                  #
                    for j from 1 do
                        iters:= iters+1:
                      #
                      # Get updated x-value
                      #
                        p:= (X__1+X__2)/2.0:
                      #
                      # Exit if required tolerance has be achieved
                      #
                        if   abs(g(p)) < err
                        then break
                      #
                      # Otherwise update the range end-points
                      #
                        elif g(p)>0
                        then X__2:=p:
                        else X__1:=p;
                        fi;
                   od:
                   return  root=p, residual=g(p), numberOfIterations=iters;
              end proc:
#
# Procedure wich computes root using Newton's method
#
  newton:= proc( g, X__A, X__B, err )
                 local j, p, X__1:=X__A, X__2:=X__B, iters:=0:
               #
               # Iterate until tolerance is achieved
               #
                 for j from 1 to 50 do
                     iters:= iters+1:
                   #
                   # Get updated x-value
                   #
                     p:= evalf(X__1-eval(g(X__1)/diff(g(x),x), x=X__1)):
                   #
                   # Exit if required tolerance has be achieved
                   #
                     if   abs(g(p)) < err
                     then break
                   #
                   # Otherwise update the
                   #
                     else X__1:=p;
                     fi;
                od:
                return  root=p, residual=g(p), numberOfIterations=iters;
          end proc:

#
# Toy example
#
   f:=x->x*exp(x)-4:
   X__1:=0: X__2:=2:
   tol:=1.0e-06:
   plot( f(x), x=0..2);
#
# Try each of these three methods
#
   regFals(f, X__1, X__2, tol);
   bisection(f, X__1, X__2, tol);
   newton(f, X__1, X__2, tol);

 

root = 1.202167739, residual = -0.983e-6, numberOfIterations = 21

 

root = 1.202167989, residual = 0.849e-6, numberOfIterations = 22

 

root = 1.202167873, residual = -0.1e-8, numberOfIterations = 9

(1)

 


Download rtFinf.mw

@AHSAN 

I have already fixed this elsewhere in this thread

@AHSAN 

and restricted the search range for the root to a region where the function returns real (rather than complex) values, produces the following

restart

with(Student[Precalculus])

k := .1

.1

(1)

F := (-(288*(lambda+2*k*(1/3)-2/3))*sqrt(2)*lambda^2*arctan((1/2)*sqrt(4*lambda/(1-k)-2)*sqrt(2))/(1-k)^2-(144*(lambda+2*k*(1/3)-2/3))*Pi*lambda^2*sqrt(2)/(1-k)^2-(36*((lambda+2*k*(1/3)-2/3)*(4*lambda/(1-k)-2)+10*lambda*(1/3)+4*k*(1/3)-4/3))*sqrt(4*lambda/(1-k)-2))*(1-k)^2/(256*lambda^2)

0.3164062500e-2*(-355.5555555*(lambda-.6000000000)*2^(1/2)*lambda^2*arctan((1/2)*(4.444444444*lambda-2)^(1/2)*2^(1/2))-558.5053607*(lambda-.6000000000)*lambda^2*2^(1/2)-36*((lambda-.6000000000)*(4.444444444*lambda-2)+(10/3)*lambda-1.200000000)*(4.444444444*lambda-2)^(1/2))/lambda^2

(2)

F := unapply(F, lambda)

proc (lambda) options operator, arrow; 0.3164062500e-2*(-355.5555555*(lambda-.6000000000)*2^(1/2)*lambda^2*arctan((1/2)*(4.444444444*lambda-2)^(1/2)*2^(1/2))-558.5053607*(lambda-.6000000000)*lambda^2*2^(1/2)-36*((lambda-.6000000000)*(4.444444444*lambda-2)+(10/3)*lambda-1.200000000)*(4.444444444*lambda-2)^(1/2))/lambda^2 end proc

(3)

seq(evalf(F(j)), j = .1 .. 1, .1)

1.249560826+2.364036895*I, .9996486610+.8792998373*I, .7497364956+.3142900773*I, .4998243302+0.5061690663e-1*I, .2056429495, -.2066756969, -.6567525583, -1.124424887, -1.601901820, -2.085401432

(4)

tol := 0.10e-5; `&lambda;__1` := .5; `&lambda;__2` := 1.0

0.10e-5

 

.5

 

1.0

(5)

for j do p := evalf([Line([`&lambda;__1`, F(`&lambda;__1`)], [`&lambda;__2`, F(`&lambda;__2`)])][4]); if abs(evalf(F(p))) < tol then break elif evalf(F(p)) > 0 then `&lambda;__2` := p else `&lambda;__1` := p end if end do; root = p, residual = evalf(F(p))

.5448797394

 

.5518826220

 

.5515854776

 

.5515867989

 

root = .5515867989, residual = -0.20e-8

(6)

``


 

Download try2.mw

@AHSAN 

reading the help at ?odetest

Upload an executable worksheet using the big green up-arrow in the Mapleprimes toolbar.

Ideally this will be sufficiently complete to illustrate the problem you are having, including defnitions for all required quantities

All Maple-related communication should take place only on this site.

You asked

I have this problem from which I want to separate constant but could not get after following your said procedure. can you pls solve for me
restart;
ode := N*diff(Phi(eta), eta, eta) - M^2*Phi(eta) = p;
                     
bc := Phi(-sigma) = 1, Phi(sigma) = -k;

The attached shows the

  1. complete solution
  2. the values of the two integration constants (although why you want these separately beats me

#
# Set up and solve the problem the easy way.
#
  ode := N*diff(Phi(eta), eta, eta) - M^2*Phi(eta) = p;
  bc := Phi(-sigma) = 1, Phi(sigma) = -k;
  sol:=simplify(dsolve([ode, bc]));

N*(diff(diff(Phi(eta), eta), eta))-M^2*Phi(eta) = p

 

Phi(-sigma) = 1, Phi(sigma) = -k

 

Phi(eta) = ((M^2+p)*exp(-M*(eta-3*sigma)/N^(1/2))+(M^2*k-p)*exp(-M*(eta-sigma)/N^(1/2))+(-M^2*k+p)*exp(M*(eta+3*sigma)/N^(1/2))+(-M^2-p)*exp(M*(eta+sigma)/N^(1/2))-p*(exp(4*M*sigma/N^(1/2))-1))/((exp(4*M*sigma/N^(1/2))-1)*M^2)

(1)

#
# Get the general solution in terms of a couple
# of integration constants
#
  genSol:=dsolve(ode);
#
# Use odetest() on this solution and the supplied
# boundary conditions, to determine required
# values for the integration constants
#
  solve
  ( odetest
    ( genSol, [ode, bc])[2..3],
    [_C1, _C2]
  )[];
#
# Substitute these constants of integration in the
# general solution - givess the saem answer as the
# previous execution group
#
  simplify( eval(genSol, %));

Phi(eta) = exp(M*eta/N^(1/2))*_C2+exp(-M*eta/N^(1/2))*_C1-p/M^2

 

[_C1 = -(M^2*exp(-M*sigma/N^(1/2))*k+M^2*exp(M*sigma/N^(1/2))+p*exp(M*sigma/N^(1/2))-exp(-M*sigma/N^(1/2))*p)/(M^2*((exp(-M*sigma/N^(1/2)))^2-(exp(M*sigma/N^(1/2)))^2)), _C2 = (k*M^2*exp(M*sigma/N^(1/2))+M^2*exp(-M*sigma/N^(1/2))-p*exp(M*sigma/N^(1/2))+exp(-M*sigma/N^(1/2))*p)/(M^2*((exp(-M*sigma/N^(1/2)))^2-(exp(M*sigma/N^(1/2)))^2))]

 

Phi(eta) = ((M^2+p)*exp(-M*(eta-3*sigma)/N^(1/2))+(M^2*k-p)*exp(-M*(eta-sigma)/N^(1/2))+(-M^2*k+p)*exp(M*(eta+3*sigma)/N^(1/2))+(-M^2-p)*exp(M*(eta+sigma)/N^(1/2))-p*(exp(4*M*sigma/N^(1/2))-1))/((exp(4*M*sigma/N^(1/2))-1)*M^2)

(2)

 

Download odeProb.mw

@AHSAN 

in the Mapleprimes toolbar to upload your complete worksheet

 

@AHSAN 

than the attached

  restart:
  with(plots):

#
# Set up a "toy" example of an 2nd order ode
# with no initial/ boundary conditions and
# solve it.
#
# Note that the colution contains two "unknown"
# constants which Maple writes as _C1 and _C2
#
  genSol:= dsolve(  diff(y(x),x,x) = 2*y(x) + 1
                 );

y(x) = exp(2^(1/2)*x)*_C2+exp(-2^(1/2)*x)*_C1-1/2

(1)

#
# So how are these two unknown constants to be
# determined? The only way is to actually provide
# two boundary conditions. So one could repeat
# the above problem, but this time with a couple
# of boundary conditions - as in the following,
# which will produce a solution with no arbitrary
# constants
#
  Sol:= dsolve( [ diff(y(x),x,x) = 2*y(x) + 1,
                   y(0)=1,
                   D(y)(0)=0
                 ]
               );

y(x) = (3/4)*exp(2^(1/2)*x)+(3/4)*exp(-2^(1/2)*x)-1/2

(2)

#
# Despite what studying textbooks might lead you,
# to believe, most ODEs cannot be solved "exactly"
# and have to be solved "numerically", which Maple
# will also do - the following will return a procedure
# from which a "numerical" solution to the original
# problem may be determined
#
  numSol:= dsolve( [ diff(y(x),x,x) = 2*y(x) + 1,
                   y(0)=1,
                   D(y)(0)=0
                 ],
                 numeric
               );

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _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](x_rkf45) else _xout := evalf(x_rkf45) 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.16825529186138485e-2, (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), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8])]), ( 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) = 3.0}, 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, X, Y, YP) option `[Y[1] = y(x), Y[2] = diff(y(x),x)]`; YP[2] := 2*Y[1]+1; YP[1] := Y[2]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(x), Y[2] = diff(y(x),x)]`; YP[2] := 2*Y[1]+1; YP[1] := Y[2]; 0 end proc, -1, 0, 0, 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 _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(0..0, {}), (3) = [x, y(x), diff(y(x), x)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["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(x_rkf45, '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(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), 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(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

(3)

#
# So now we can plot the "exact" solution and the
# "numerical" solution on the saem graph, just to
# see if the difference is noticeable - and for
# this ODE, it really isn't
#
  display
  ( [ plot
      ( rhs(Sol),
        x=0..1
      ),
      odeplot
      ( numSol,
        [x,y(x)],
        x=0..1
      )
    ],
    color=[red, blue],
    title="Exact and Numerical Solutions"
  );
#
# Although if one really wants to try one could
# actually plot the difference between these two
# solutions, just to see how accurate the numerical
# solution process actually is. Before you get excited
# by this plot, note the vertical scale indicating
# that the error between the exact and numerical
# solutions is O(10^-7)
#
  odeplot
  ( numSol,
    [x, y(x)-rhs(Sol)],
    x=0..1,
    title="Error between exact and numeric solutions"
  );

 

 

#
# An alternative method of determining the values of the
# constants in the general solution - ie genSol above, is
# to use the odetest() command. Essentially this substitutes
# the general solution, with the unknowns (_C1 and _C2)
# into the original oded and supply the same boundary
# conditions. Since all entries in the returned list should
# be zero, the second and third entries in the returned list
# give a pair of equations whihc the unknown constants _C1
# and _C2 must satisfy
#
  otest:= odetest( genSol,
                   [ diff(y(x),x,x) = 2*y(x) + 1,
                     y(0)=1,
                     D(y)(0)=0
                   ]
                 );

[0, 3/2-_C2-_C1, -2^(1/2)*_C2+2^(1/2)*_C1]

(4)

#
# Determine the values of the integration constants
# which make all above list entries zero
#
  cVals:=solve(otest[-2..-1]);

{_C1 = 3/4, _C2 = 3/4}

(5)

#
# Now substitute the values of these constants in the
# general solution. Note that this is identical to
# the exact solution 'Sol' obtained above
#
  subs( cVals, genSol);

y(x) = (3/4)*exp(2^(1/2)*x)+(3/4)*exp(-2^(1/2)*x)-1/2

(6)

 


 

Download solvingODEs.mw

Maple 2020 returns the two solutions you (and dharr) have already obtained, with no warniings. I wouldn't like to guarantee that this means there are no other solutions - probably just that the existence of any other solutions is very very unlikely.

restart;
interface(version);
solve({l*(2*l^2*lambda^4*sigma*w*a[2]+l^2*lambda^2*mu*w*b[1]+6*l*lambda^2*m*sigma*a[0]^2-6*l*lambda^2*m*b[1]^2+6*l*m*mu^2*a[0]^2-l*lambda^2*    rho*sigma*a[0]-l*mu^2*rho*a[0]+4*lambda^2*sigma*w*a[0]+4*mu^2*w*a[0]) = 0,
 l*(2*l^2*lambda^3*sigma*w*a[1]+6*l^2*lambda^2*mu*w*b[2]+2*l^2*lambda*mu^2*w*a[1]+12*l*lambda^2*m*sigma*a[0]*a[1]-12*l*lambda^2*m*b[1]*b[2]-l*lambda^2*rho*sigma*a[1]+12*l*m*mu^2*a[0]*a[1]-l*mu^2*rho*a[1]+4*lambda^2*sigma*w*a[1]+4*mu^2*w*a[1]) = 0,
 l*(5*l^2*lambda^3*sigma*w*b[2]-3*l^2*lambda^2*mu*sigma*w*a[1]-7*l^2*lambda*mu^2*w*b[2]-3*l^2*mu^3*w*a[1]+12*l*lambda^2*m*sigma*a[0]*b[2]+12*l*lambda^2*m*sigma*a[1]*b[1]-l*lambda^2*rho*sigma*b[2]+24*l*lambda*m*mu*b[1]*b[2]+12*l*m*mu^2*a[0]*b[2]+12*l*m*mu^2*a[1]*b[1]-l*mu^2*rho*b[2]+4*lambda^2*sigma*w*b[2]+4*mu^2*w*b[2]) = 0,
 l*(8*l^2*lambda^3*sigma*w*a[2]+6*l^2*lambda*mu^2*w*a[2]+12*l*lambda^2*m*sigma*a[0]*a[2]+6*l*lambda^2*m*sigma*a[1]^2+l^2*lambda*mu*w*b[1]-6*l*lambda^2*m*b[2]^2-l*lambda^2*rho*sigma*a[2]+12*l*m*mu^2*a[0]*a[2]+6*l*m*mu^2*a[1]^2-6*l*lambda*m*b[1]^2-l*mu^2*rho*a[2]+4*lambda^2*sigma*w*a[2]+4*mu^2*w*a[2]) = 0,
 -l*(4*l^2*lambda^3*mu*sigma*w*a[2]-l^2*lambda^3*sigma*w*b[1]+l^2*lambda*mu^2*w*b[1]-12*l*lambda^2*m*sigma*a[0]*b[1]+l*lambda^2*rho*sigma*b[1]-12*l*lambda*m*mu*b[1]^2-12*l*m*mu^2*a[0]*b[1]+l*mu^2*rho*b[1]-4*lambda^2*sigma*w*b[1]-4*mu^2*w*b[1]) = 0,
 6*l^2*(l*lambda^2*sigma*w*a[2]+lambda^2*m*sigma*a[2]^2+l*mu^2*w*a[2]+m*mu^2*a[2]^2-lambda*m*b[2]^2) = 0, 2*l^2*(l*lambda^2*sigma*w*a[1]+6*lambda^2*m*sigma*a[1]*a[2]+3*l*lambda*mu*w*b[2]+l*mu^2*w*a[1]+6*m*mu^2*a[1]*a[2]-6*lambda*m*b[1]*b[2]) = 0,
 -2*l^2*(5*l*lambda^2*mu*sigma*w*a[2]-l*lambda^2*sigma*w*b[1]+5*l*mu^3*w*a[2]-6*lambda^2*m*sigma*a[1]*b[2]-6*lambda^2*m*sigma*a[2]*b[1]-l*mu^2*w*b[1]-6*lambda*m*mu*b[2]^2-6*m*mu^2*a[1]*b[2]-6*m*mu^2*a[2]*b[1]) = 0,
 6*l^2*b[2]*(l*w+2*m*a[2]) = 0}, {a[0], a[1], a[2], b[1], b[2]});

`Standard Worksheet Interface, Maple 2020.0, Windows 7, March 4 2020 Build ID 1455132`

 

{a[0] = 0, a[1] = 0, a[2] = 0, b[1] = 0, b[2] = 0}, {a[0] = (1/6)*(l*rho-4*w)/(l*m), a[1] = 0, a[2] = 0, b[1] = 0, b[2] = 0}

(1)

 

Download polySol.mw

 

@ogunmiloro 

on your machine you will have to change the line

someData:=ExcelTools:-Import("C:/Users/TomLeslie/Desktop/testDAta.xls"):

to something that makes sense for your machine.

The filePath in the above command, ie

"C:/Users/TomLeslie/Desktop/testDAta.xls"

only refers to the location + filename I used on my computer to store the data for test purposes.

I have no idea where the data file is stored on your machine! Only you know this!

@raj2018 

If you examine the output of the dsolve() command in my original post, you will note that it contains q1(X) = q1(X) - in other words q1(x) is entirely arbitrary.

You can pick any function you want for q1(x), substitute it in the expressions for all the other functions, and the solution is still valid

@HaHu 

The concept of recursion is very useful, but only when there is a "terminating case". So

x:=10
x:=x+1;

is perfectly fine, because the right-hand-side 'x+1', will evaluate to '10+1', ie 11, and 11 is assigned to 'x'. No problem

However

p:=p+1

doesn't make much sense. If applied in a brain dead way the right hand side 'p+1' will become '(p+1)+1' which will evaluate to '((p+1)+1)+1) - and so on ad infinitum

Asssignments are not equations - it is obviously true that one could write the equation  p=p2+1 and "solve" this equation for p, but the assignment p:=p2+1 doesn't make any sense

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