Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 166 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

More of an observation really...

tomleslie 13876
October 04 2020
1 3

Looks like some "order-of-calculation" issue. If you supply 'deps' as a set rather than a list, then both dsolve(sys) and dsolve(sys,deps) give the same answer. Check the attached


 

> 

restart;
sys:=[ diff(x(t),t) = 2*x(t)-z(t),
       diff(y(t),t) = 2*y(t)+z(t),
       diff(z(t),t)=2*z(t),
       diff(w(t),t)=-z(t)+2*w(t)
     ];
deps:={x(t),y(t),z(t),w(t)};

dsolve(sys);

dsolve(sys,deps);

[diff(x(t), t) = 2*x(t)-z(t), diff(y(t), t) = 2*y(t)+z(t), diff(z(t), t) = 2*z(t), diff(w(t), t) = -z(t)+2*w(t)]

  

{w(t), x(t), y(t), z(t)}

  

{w(t) = (-_C4*t+_C3)*exp(2*t), x(t) = (-_C4*t+_C2)*exp(2*t), y(t) = (_C4*t+_C1)*exp(2*t), z(t) = _C4*exp(2*t)}

  

{w(t) = (-_C4*t+_C3)*exp(2*t), x(t) = (-_C4*t+_C2)*exp(2*t), y(t) = (_C4*t+_C1)*exp(2*t), z(t) = _C4*exp(2*t)}

 (1)
> 

 


 

Download listset.mw

Maybe...

tomleslie 13876
October 03 2020
1 0

There are (at least) three possible methds to solve this problem

  1. Direct integration as suggested by Preben - which is defintely the simplest
  2. Using dsolve() with no boundary conditions to return an analytic solution containing an integration constant. Take the limit of this solution as x->-infinity, to determine the value of the integration constant which enures that this limit is 0
  3. Using a zero-finding process. Basically,
    1. set the initial condition P(0)=c:
    2. solve the ode numerically and obtain P(-10^9) as a proxy for P(-infinity)
    3. use fsolve() to determine the value of c, such that P(-10^9)=0
    4. having obtained the appropriate value of 'c', solve the ODE numerically with the boundary condition P(0)=c
    5. This process is a "shooting method" (more or less)

The attached shows methods (2) and (3) in the above

> 

  restart;
> 

  sigma := x^2 + 1;
  M := 1/10;
  N := 1/2;
  ode := diff(P(x), x) = -3*(M - 1)/(2*sigma^2) - 3*N/sigma^3;

x^2+1

  

1/10

  

1/2

  

diff(P(x), x) = (27/20)/(x^2+1)^2-(3/2)/(x^2+1)^3

 (1)
> 

#
# Solve the ODE analytically, with no boundary conditions
# (Solution will include a constant)
#
  sol := dsolve(ode):
#
# Determine the value of the  constant so that P(-infinity)=0
# and substitute this in the solution
#
  sol := subs(isolate(limit(rhs(sol), x = -infinity) = 0, _C1), sol);
#
# Plot this solution
#
  plot(rhs(sol), x = -10 .. 10);

P(x) = (3/20)*((3/4)*x^3-(7/4)*x)/(x^2+1)^2+(9/80)*arctan(x)+(9/160)*Pi

  

 
> 

#
# Procedure which accepts a parameter 'c' which is used
# in a boundary condition P(0)=c, and returns the value
# P(-10^9) which is a proxy for P(-infinity)
#
  S:= proc(c)
          return rhs
                 ( dsolve
                   ( [ode, P(0)=c],
                     numeric
                   )(-10^9)[2]
                 );
      end proc:
#
# Use fsolve to determine the value of the parameter c
# in the condition P(0)=c, such that the condition
# P(-10^9)=0 is obtained
#
  ic:=fsolve( 'S'(a)=0);
#
# Numerically solve the ode with the boundary condition
# P(0)=ic, where the value of ic (determined above)
# ensures that P(-10^9)=0
#
  sol:=dsolve( [ode, P(0)=ic], numeric);
#
# Plot this numeric solution to compare with the analytic
# solution obtained previously
#
  plots:-odeplot( sol, [x, P(x)], x=-10..10);

.1767146001

  

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) = 1, (2) = 1, (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.9311739160426387e-1, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = .1767146001}, 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..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = .1767146001}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -.1499999999999999}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(x)]`; YP[1] := (27/20)/(X^2+1)^2-(3/2)/(X^2+1)^3; 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] = P(x)]`; YP[1] := (27/20)/(X^2+1)^2-(3/2)/(X^2+1)^3; 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..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _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, P(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

  

 
> 

 


 

Download shoot.mw

Hmmm...

tomleslie 13876
October 02 2020
0 3

As a "general-puprose" implementation of Newton's method I can see all sorts of potential issues with this code - but if I just "make it work", with minimal changes, see the attached

> 

restart;
Newton:=proc(f, a)
             local n,
                   x1:= a-f(a)/eval( diff(f(x),x),x=a);
             for n to 9 do
                 x1:=x1-f(x1)/eval(diff(f(x),x), x=x1);
             end do;
             if abs( f(x1) ) <= 10^(-8)
             then return x1;
             else return FAIL
             fi;
        end proc:
f:=x->x^3-3*x^2+2*exp(x)+5;
Newton( f, -1.2);

proc (x) options operator, arrow; x^3-3*x^2+2*exp(x)+5 end proc

  

-1.162515305

 (1)
> 

#
# Check the above
#
  fsolve(f(x));

-1.162515305

 (2)
> 

 

Download newt.mwDownload newt.mw

Use the FileTools package...

tomleslie 13876
October 02 2020
1 0

I agree it looks like a bug, but the attached does work without any string conversions (Change the file path to somethng appropriate for your machine)

a := 2^1000

10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

 (1)

FileTools[Text]:-WriteInteger("C:/Users/TomLeslie/Desktop/a.txt", a); FileTools[Text]:-Close("C:/Users/TomLeslie/Desktop/a.txt")

restart

FileTools[Text]:-ReadInteger("C:/Users/TomLeslie/Desktop/a.txt")

10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

 (2)

NULL


 

Download readInt.mw

plot rhs(sol)...

tomleslie 13876
September 30 2020
3 1

as in the atached

> 

restart;
pde := diff(u(x, t), t $ 2) - diff(u(x, t), x $ 2) = 0;
h := piecewise(-1 <= x and x <= 0, x + 2, 0 <= x and x <= 1, 2 - x, 0);
ic := u(x, 0) = h, D[2](u)(x, 0) = 0;
sol := pdsolve([pde, ic]);
plots[animate](plot, [rhs(sol), x = -10 .. 10], t = 0 .. 5)

pde := diff(u(x, t), t, t)-(diff(u(x, t), x, x)) = 0

  

h := piecewise(-1 <= x and x <= 0, x+2, 0 <= x and x <= 1, 2-x, 0)

  

ic := u(x, 0) = piecewise(-1 <= x and x <= 0, x+2, 0 <= x and x <= 1, 2-x, 0), (D[2](u))(x, 0) = 0

  

u(x, t) = (1/2)*piecewise(-t+x < -1, 0, -t+x <= 0, -t+x+2, -t+x <= 1, 2+t-x, 1 < -t+x, 0)+(1/2)*piecewise(t+x < -1, 0, t+x <= 0, t+x+2, t+x <= 1, 2-t-x, 1 < t+x, 0)

  

 
> 

 


 

Download anim.mw

Well...

tomleslie 13876
September 29 2020
2 1

You must not use [] parentheses to group terms in an expression - these are used to signify indexable quantities (lists, matrices etc)

Your expression contains the variable 'y' which is nowhere defined.

You should post code rather than pictures of code because most people here won't retype your code from scratch (and the process is error-prone so I'm not guaranteeing that I retyped correctly!!). Use the big green up-arrow in the Mapleprimes toolbar to upload your worksheet

If I correct the '[]' issue and arbitrarily set y=1, I get the attached

> 

restart;
> 

with(plots):
with(DETools):
b:=0.2:
m:=0.1:
y:=1:
> 

DE:=diff(x(N),N)=1/2*x(N)*(2*(3*(b*(x(N)^2)^2-4*m^2+x(N)^2+4))/(2*m^2+y^2-2)+3*(b*x(N))^2+5);

diff(x(N), N) = (1/2)*x(N)*(-1.224489796*x(N)^4-19.24489796-6.002448978*x(N)^2)

 (1)
> 

DEplot(DE, x(N), N=-15..5, x=0..1);

 
> 

 


 

Download dplt.mw

 

Definitely interesting...

tomleslie 13876
September 29 2020
0 0

because performing the calsulation in two slightly different ways gives two somewhat different outcomes. See the comments in the attached

> 

restart;
> 

#
# With this definition, the integral is not solved
# until a numerical value for the passed parameter
# is supplied. This numerical value is substituted
# into the integrand, and the value of the integral
# is then computed
#
  expr1:=a->int(exp(-(x^2-a)^2),x=0..infinity);
  plot(expr1(a), a=-2..5);
#
# Check some values, including a=0
#
  seq(evalf(expr1(j)), j in [-1, -1/2, 0, 1/2, 1]);

proc (a) options operator, arrow; int(exp(-(x^2-a)^2), x = 0 .. infinity) end proc

  

 

.2059311656, .5328659528, .9064024772, 1.075270372, .9868660750

 (1)
> 

#
# With this method, the integral is solved for symbolic
# values of the parameter 'a'. All the subsequent
# evaluations susbtitute 'a' into the solution of the
# integral
#
  expr2:=unapply(int(exp(-(x^2-a)^2),x=0..infinity),a);
  plot(expr2(a), a=-2..5);
#
# Check some values. Note that for an argument of '0'
# evaluation 'fails', (presumably because of the
# existence of 'a' in the denominator)
#
  seq(evalf(expr2(j)), j in [-1, -1/2, 1/2, 1]);
  expr2(0);
#
# So try evaluating at zero, b using a limit
#
  evalf(limit( expr2(x), x=0));

proc (a) options operator, arrow; (1/4)*exp(-a^2)*(a^2)^(1/4)*2^(1/2)*exp((1/2)*a^2)*(BesselK(7/4, (1/2)*a^2)*a^2-3*BesselK(3/4, (1/2)*a^2))/a^2 end proc

  

 

.2059311656, .5328659615, .5328659615, .2059311656

  

Error, (in BesselK) numeric exception: division by zero

  

.9064024772

 (2)
> 

 

> 

 


 

Download badInt.mw

Try the 'approxsoln' option...

tomleslie 13876
September 24 2020
0 4

The attached solves the system for lambda=0.1. It then uses this solution (via approxsoln=..) to obtain a solution for lambsa=0.2


 

restart

b := 25.4

f := 2.2

h := 10

E := 0.70e5

G := 0.26e5

Ec := 420

Gc := 220

nuc := .25

A11 := E*f

D11 := (1/12)*E*f^3

A55 := G*f

dt := fNULL

db := f

c := h

phi_ct := .55; delta0_ct := .12; phi_cb := .55; delta0_cb := .12

eq1_1 := (diff(y1_1(x), x))/(27.1) = 0

eq1_2 := (diff(y1_2(x), x))/(27.1) = 0

eq1_3 := (diff(y1_3(x), x))/(27.1) = y1_2(x)

eq1_4 := (diff(y1_4(x), x))/(27.1) = 2*b*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2

eq1_5 := (diff(y1_5(x), x))/(27.1) = b*phi_cb*(y1_18(x)-y1_24(x)-y1_25(x)*c-y1_26(x)*c^2)/delta0_cb^2

eq1_6 := (diff(y1_6(x), x))/(27.1) = y1_5(x)-2*b*db*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2*(1/2)

eq1_7 := (diff(y1_7(x), x))/(27.1) = -2*b*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2

eq1_8 := (diff(y1_8(x), x))/(27.1) = -2*b*c*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2+y1_11(x)

eq1_9 := (diff(y1_9(x), x))/(27.1) = -2*b*c^2*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2+2*y1_12(x)

eq1_10 := (diff(y1_10(x), x))/(27.1) = -2*b*c^3*phi_cb*(y1_17(x)+(1/2)*db*y1_19(x)-y1_20(x)-y1_21(x)*c-y1_22(x)*c^2-y1_23(x)*c^3)/delta0_cb^2+3*y1_13(x)

eq1_11 := (diff(y1_11(x), x))/(27.1) = -b*phi_cb*(y1_18(x)-y1_24(x)-y1_25(x)*c-y1_26(x)*c^2)/delta0_cb^2

eq1_12 := (diff(y1_12(x), x))/(27.1) = -b*c*phi_cb*(y1_18(x)-y1_24(x)-y1_25(x)*c-y1_26(x)*c^2)/delta0_cb^2+2*Ec*b*c*y1_25(x)+nuc*y1_7(x)

eq1_13 := (diff(y1_13(x), x))/(27.1) = -b*c^2*phi_cb*(y1_18(x)-y1_24(x)-y1_25(x)*c-y1_26(x)*c^2)/delta0_cb^2+2*(4*Ec*b*c^3*y1_26(x)*(1/3)+nuc*y1_8(x))

eq1_14 := (diff(y1_14(x), x))/(27.1) = y1_1(x)/(b*A11)

eq1_15 := (diff(y1_15(x), x))/(27.1) = y1_16(x)+y1_2(x)/(b*A55)

eq1_16 := (diff(y1_16(x), x))/(27.1) = -y1_3(x)/(b*D11)

eq1_17 := (diff(y1_17(x), x))/(27.1) = y1_4(x)/(b*A11) 

eq1_18 := (diff(y1_18(x), x))/(27.1) = y1_19(x)+y1_5(x)/(b*A55)

eq1_19 := (diff(y1_19(x), x))/(27.1) = -y1_6(x)/(b*D11)

eq1_20 := (diff(y1_20(x), x))/(27.1) = (-nuc^2+1)*(9*c^2*y1_7(x)-15*y1_9(x))/(8*Ec*b*c^3)-nuc*y1_25(x)

eq1_21 := (diff(y1_21(x), x))/(27.1) = (-nuc^2+1)*(75*c^2*y1_8(x)-105*y1_10(x))/(8*Ec*b*c^5)-2*nuc*y1_26(x)

eq1_22 := (diff(y1_22(x), x))/(27.1) = (-nuc^2+1)*(45*y1_9(x)-15*c^2*y1_7(x))/(8*Ec*b*c^5)

eq1_23 := (diff(y1_23(x), x))/(27.1) = (-nuc^2+1)*(175*y1_10(x)-105*c^2*y1_8(x))/(8*Ec*b*c^7)

eq1_24 := (diff(y1_24(x), x))/(27.1) = (9*c^2*y1_11(x)-15*y1_13(x))/(8*Gc*b*c^3)-y1_21(x)

eq1_25 := (diff(y1_25(x), x))/(27.1) = 3*y1_12(x)/(2*Gc*b*c^3)-2*y1_22(x)

eq1_26 := (diff(y1_26(x), x))/(27.1) = (45*y1_13(x)-15*c^2*y1_11(x))/(8*Gc*b*c^5)-3*y1_23(x)

NULL

eq2_1 := (diff(y2_1(x), x))/(23.7) = 0; eq2_2 := (diff(y2_2(x), x))/(23.7) = 0; eq2_3 := (diff(y2_3(x), x))/(23.7) = y2_2(x); eq2_4 := (diff(y2_4(x), x))/(23.7) = 2*b*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2; eq2_5 := (diff(y2_5(x), x))/(23.7) = b*phi_cb*(y2_18(x)-y2_24(x)-y2_25(x)*c-y2_26(x)*c^2)/delta0_cb^2; eq2_6 := (diff(y2_6(x), x))/(23.7) = y2_5(x)-2*b*db*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2*(1/2); eq2_7 := (diff(y2_7(x), x))/(23.7) = -2*b*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2; eq2_8 := (diff(y2_8(x), x))/(23.7) = -2*b*c*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2+y2_11(x); eq2_9 := (diff(y2_9(x), x))/(23.7) = -2*b*c^2*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2+2*y2_12(x); eq2_10 := (diff(y2_10(x), x))/(23.7) = -2*b*c^3*phi_cb*(y2_17(x)+(1/2)*db*y2_19(x)-y2_20(x)-y2_21(x)*c-y2_22(x)*c^2-y2_23(x)*c^3)/delta0_cb^2+3*y2_13(x); eq2_11 := (diff(y2_11(x), x))/(23.7) = -b*phi_cb*(y2_18(x)-y2_24(x)-y2_25(x)*c-y2_26(x)*c^2)/delta0_cb^2; eq2_12 := (diff(y2_12(x), x))/(23.7) = -b*c*phi_cb*(y2_18(x)-y2_24(x)-y2_25(x)*c-y2_26(x)*c^2)/delta0_cb^2+2*Ec*b*c*y2_25(x)+nuc*y2_7(x); eq2_13 := (diff(y2_13(x), x))/(23.7) = -b*c^2*phi_cb*(y2_18(x)-y2_24(x)-y2_25(x)*c-y2_26(x)*c^2)/delta0_cb^2+2*(4*Ec*b*c^3*y2_26(x)*(1/3)+nuc*y2_8(x)); eq2_14 := (diff(y2_14(x), x))/(23.7) = y2_1(x)/(b*A11); eq2_15 := (diff(y2_15(x), x))/(23.7) = y2_16(x)+y2_2(x)/(b*A55); eq2_16 := (diff(y2_16(x), x))/(23.7) = -y2_3(x)/(b*D11); eq2_17 := (diff(y2_17(x), x))/(23.7) = y2_4(x)/(b*A11); eq2_18 := (diff(y2_18(x), x))/(23.7) = y2_19(x)+y2_5(x)/(b*A55); eq2_19 := (diff(y2_19(x), x))/(23.7) = -y2_6(x)/(b*D11); eq2_20 := (diff(y2_20(x), x))/(23.7) = (-nuc^2+1)*(9*c^2*y2_7(x)-15*y2_9(x))/(8*Ec*b*c^3)-nuc*y2_25(x); eq2_21 := (diff(y2_21(x), x))/(23.7) = (-nuc^2+1)*(75*c^2*y2_8(x)-105*y2_10(x))/(8*Ec*b*c^5)-2*nuc*y2_26(x); eq2_22 := (diff(y2_22(x), x))/(23.7) = (-nuc^2+1)*(45*y2_9(x)-15*c^2*y2_7(x))/(8*Ec*b*c^5); eq2_23 := (diff(y2_23(x), x))/(23.7) = (-nuc^2+1)*(175*y2_10(x)-105*c^2*y2_8(x))/(8*Ec*b*c^7); eq2_24 := (diff(y2_24(x), x))/(23.7) = (9*c^2*y2_11(x)-15*y2_13(x))/(8*Gc*b*c^3)-y2_21(x); eq2_25 := (diff(y2_25(x), x))/(23.7) = 3*y2_12(x)/(2*Gc*b*c^3)-2*y2_22(x); eq2_26 := (diff(y2_26(x), x))/(23.7) = (45*y2_13(x)-15*c^2*y2_11(x))/(8*Gc*b*c^5)-3*y2_23(x)

NULL

eq3_1 := (diff(y3_1(x), x))/(101.6) = -2*b*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2

eq3_2 := (diff(y3_2(x), x))/(101.6) = -b*phi_ct*(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2

eq3_3 := (diff(y3_3(x), x))/(101.6) = y3_2(x)-2*b*dt*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2*(1/2)

eq3_4 := (diff(y3_4(x), x))/(101.6) = 2*b*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2

eq3_5 := (diff(y3_5(x), x))/(101.6) = b*phi_cb*(y3_18(x)-y3_24(x)-y3_25(x)*c-y3_26(x)*c^2)/delta0_cb^2

eq3_6 := (diff(y3_6(x), x))/(101.6) = y3_5(x)-2*b*db*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2*(1/2)

eq3_7 := (diff(y3_7(x), x))/(101.6) = -2*b*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2+2*b*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2

eq3_8 := (diff(y3_8(x), x))/(101.6) = -2*b*c*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2-2*b*c*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2+y3_11(x)

eq3_9 := (diff(y3_9(x), x))/(101.6) = -2*b*c^2*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2+2*b*c^2*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2+2*y3_12(x)

eq3_10 := (diff(y3_10(x), x))/(101.6) = -2*b*c^3*phi_ct*(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))*(1+(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct)*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2-2*b*c^3*phi_cb*(y3_17(x)+(1/2)*db*y3_19(x)-y3_20(x)-y3_21(x)*c-y3_22(x)*c^2-y3_23(x)*c^3)/delta0_cb^2+3*y3_13(x)

eq3_11 := (diff(y3_11(x), x))/(101.6) = -b*phi_cb*(y3_18(x)-y3_24(x)-y3_25(x)*c-y3_26(x)*c^2)/delta0_cb^2+b*phi_ct*(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2

eq3_12 := (diff(y3_12(x), x))/(101.6) = -b*c*phi_ct*(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2-b*c*phi_cb*(y3_18(x)-y3_24(x)-y3_25(x)*c-y3_26(x)*c^2)/delta0_cb^2+2*Ec*b*c*y3_25(x)+nuc*y3_7(x)

eq3_13 := (diff(y3_13(x), x))/(101.6) = -b*c^2*phi_cb*(y3_18(x)-y3_24(x)-y3_25(x)*c-y3_26(x)*c^2)/delta0_cb^2+b*c^2*phi_ct*(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))*exp(-(y3_24(x)-y3_25(x)*c+y3_26(x)*c^2-y3_15(x))/delta0_ct-(y3_20(x)-y3_21(x)*c+y3_22(x)*c^2-y3_23(x)*c^3-y3_14(x)+(1/2)*dt*y3_16(x))^2/delta0_ct^2)/delta0_ct^2+2*(4*Ec*b*c^3*y3_26(x)*(1/3)+nuc*y3_8(x))

eq3_14 := (diff(y3_14(x), x))/(101.6) = y3_1(x)/(b*A11); eq3_15 := (diff(y3_15(x), x))/(101.6) = y3_16(x)+y3_2(x)/(b*A55); eq3_16 := (diff(y3_16(x), x))/(101.6) = -y3_3(x)/(b*D11); eq3_17 := (diff(y3_17(x), x))/(101.6) = y3_4(x)/(b*A11); eq3_18 := (diff(y3_18(x), x))/(101.6) = y3_19(x)+y3_5(x)/(b*A55); eq3_19 := (diff(y3_19(x), x))/(101.6) = -y3_6(x)/(b*D11); eq3_20 := (diff(y3_20(x), x))/(101.6) = (-nuc^2+1)*(9*c^2*y3_7(x)-15*y3_9(x))/(8*Ec*b*c^3)-nuc*y3_25(x); eq3_21 := (diff(y3_21(x), x))/(101.6) = (-nuc^2+1)*(75*c^2*y3_8(x)-105*y3_10(x))/(8*Ec*b*c^5)-2*nuc*y3_26(x); eq3_22 := (diff(y3_22(x), x))/(101.6) = (-nuc^2+1)*(45*y3_9(x)-15*c^2*y3_7(x))/(8*Ec*b*c^5); eq3_23 := (diff(y3_23(x), x))/(101.6) = (-nuc^2+1)*(175*y3_10(x)-105*c^2*y3_8(x))/(8*Ec*b*c^7); eq3_24 := (diff(y3_24(x), x))/(101.6) = (9*c^2*y3_11(x)-15*y3_13(x))/(8*Gc*b*c^3)-y3_21(x); eq3_25 := (diff(y3_25(x), x))/(101.6) = 3*y3_12(x)/(2*Gc*b*c^3)-2*y3_22(x); eq3_26 := (diff(y3_26(x), x))/(101.6) = (45*y3_13(x)-15*c^2*y3_11(x))/(8*Gc*b*c^5)-3*y3_23(x)

NULL

eq1_27 := diff(y1_27(x), x) = 0

``

NULL

odesys := [eq1_1, eq1_2, eq1_3, eq1_4, eq1_5, eq1_6, eq1_7, eq1_8, eq1_9, eq1_10, eq1_11, eq1_12, eq1_13, eq1_14, eq1_15, eq1_16, eq1_17, eq1_18, eq1_19, eq1_20, eq1_21, eq1_22, eq1_23, eq1_24, eq1_25, eq1_26, eq1_27, eq2_1, eq2_2, eq2_3, eq2_4, eq2_5, eq2_6, eq2_7, eq2_8, eq2_9, eq2_10, eq2_11, eq2_12, eq2_13, eq2_14, eq2_15, eq2_16, eq2_17, eq2_18, eq2_19, eq2_20, eq2_21, eq2_22, eq2_23, eq2_24, eq2_25, eq2_26, eq3_1, eq3_2, eq3_3, eq3_4, eq3_5, eq3_6, eq3_7, eq3_8, eq3_9, eq3_10, eq3_11, eq3_12, eq3_13, eq3_14, eq3_15, eq3_16, eq3_17, eq3_18, eq3_19, eq3_20, eq3_21, eq3_22, eq3_23, eq3_24, eq3_25, eq3_26, y1_1(0) = 0, y1_2(0) = 0, y1_3(0) = 0, y1_4(0) = 0, y1_5(0) = 0, y1_6(0) = 0, y1_7(0) = 0, y1_8(0) = 0, y1_9(0) = 0, y1_10(0) = 0, y1_11(0) = 0, y1_12(0) = 0, y1_13(0) = 0, y1_2(1)-y2_2(0)-y1_27 = 0, y1_3(1)-y2_3(0) = 0, y1_6(1)-y2_6(0) = 0, y1_7(1)-y2_7(0) = 0, y1_8(1)-y2_8(0) = 0, y1_9(1)-y2_9(0) = 0, y1_10(1)-y2_10(0) = 0, y1_11(1)-y2_11(0) = 0, y1_12(1)-y2_12(0) = 0, y1_13(1)-y2_13(0) = 0, y1_14(1) = y2_14(0), y1_15(1) = y2_15(0), y1_16(1) = y2_16(0), y1_17(1) = y2_17(0), y1_18(1) = y2_18(0), y1_19(1) = y2_19(0), y1_20(1) = y2_20(0), y1_21(1) = y2_21(0), y1_22(1) = y2_22(0), y1_23(1) = y2_23(0), y1_24(1) = y2_24(0), y1_25(1) = y2_25(0), y1_26(1) = y2_26(0), y1_14(1) = 0, y1_15(1) = -lambda, y1_17(1) = 0, y1_18(1) = 0, y2_1(1)-y3_1(0) = 0, y2_2(1)-y3_2(0) = 0, y2_3(1)-y3_3(0) = 0, y2_4(1)-y3_4(0) = 0, y2_5(1)-y3_5(0) = 0, y2_6(1)-y3_6(0) = 0, y2_7(1)-y3_7(0) = 0, y2_8(1)-y3_8(0) = 0, y2_9(1)-y3_9(0) = 0, y2_10(1)-y3_10(0) = 0, y2_11(1)-y3_11(0) = 0, y2_12(1)-y3_12(0) = 0, y2_13(1)-y3_13(0) = 0, y2_14(1) = y3_14(0), y2_15(1) = y3_15(0), y2_16(1) = y3_16(0), y2_17(1) = y3_17(0), y2_18(1) = y3_18(0), y2_19(1) = y3_19(0), y2_20(1) = y3_20(0), y2_21(1) = y3_21(0), y2_22(1) = y3_22(0), y2_23(1) = y3_23(0), y2_24(1) = y3_24(0), y2_25(1) = y3_25(0), y2_26(1) = y3_26(0), y3_1(1) = 0, y3_2(1) = 0, y3_3(1) = 0, y3_4(1) = 0, y3_5(1) = 0, y3_6(1) = 0, y3_7(1) = 0, y3_8(1) = 0, y3_9(1) = 0, y3_10(1) = 0, y3_11(1) = 0, y3_12(1) = 0, y3_13(1) = 0], [y1_1(x), y1_2(x), y1_3(x), y1_4(x), y1_5(x), y1_6(x), y1_7(x), y1_8(x), y1_9(x), y1_10(x), y1_11(x), y1_12(x), y1_13(x), y1_14(x), y1_15(x), y1_16(x), y1_17(x), y1_18(x), y1_19(x), y1_20(x), y1_21(x), y1_22(x), y1_23(x), y1_24(x), y1_25(x), y1_26(x), y1_27(x), y2_1(x), y2_2(x), y2_3(x), y2_4(x), y2_5(x), y2_6(x), y2_7(x), y2_8(x), y2_9(x), y2_10(x), y2_11(x), y2_12(x), y2_13(x), y2_14(x), y2_15(x), y2_16(x), y2_17(x), y2_18(x), y2_19(x), y2_20(x), y2_21(x), y2_22(x), y2_23(x), y2_24(x), y2_25(x), y2_26(x), y3_1(x), y3_2(x), y3_3(x), y3_4(x), y3_5(x), y3_6(x), y3_7(x), y3_8(x), y3_9(x), y3_10(x), y3_11(x), y3_12(x), y3_13(x), y3_14(x), y3_15(x), y3_16(x), y3_17(x), y3_18(x), y3_19(x), y3_20(x), y3_21(x), y3_22(x), y3_23(x), y3_24(x), y3_25(x), y3_26(x)]
``

solnumeric1 := dsolve(eval(odesys, lambda = .1), numeric)

solnumeric1(1)

[x = 1., y1_1(x) = HFloat(0.0), y1_2(x) = HFloat(-3.420984483671309e-31), y1_3(x) = HFloat(-9.266873428961157e-30), y1_4(x) = HFloat(-2.8265796778082293), y1_5(x) = HFloat(4.585656258526918), y1_6(x) = HFloat(29.73356557713309), y1_7(x) = HFloat(2.8265796778082293), y1_8(x) = HFloat(1.6414688476893966), y1_9(x) = HFloat(78.13532261272604), y1_10(x) = HFloat(134.29458529790443), y1_11(x) = HFloat(-4.585656258526918), y1_12(x) = HFloat(-12.84587679453849), y1_13(x) = HFloat(-141.040214445319), y1_14(x) = HFloat(0.0), y1_15(x) = HFloat(-0.09999999999999995), y1_16(x) = HFloat(0.004086375487434791), y1_17(x) = HFloat(0.0), y1_18(x) = HFloat(0.0), y1_19(x) = HFloat(-2.458791784171863e-4), y1_20(x) = HFloat(-0.001957542333257033), y1_21(x) = HFloat(1.9956662771272542e-4), y1_22(x) = HFloat(-1.5079719874448682e-6), y1_23(x) = HFloat(3.3273624675095036e-8), y1_24(x) = HFloat(-7.164450668824305e-4), y1_25(x) = HFloat(1.1943543119401487e-5), y1_26(x) = HFloat(1.0902056689470582e-6), y1_27(x) = HFloat(-11.236234112909697), y2_1(x) = HFloat(7.740295426479526e-11), y2_2(x) = HFloat(11.236234112909697), y2_3(x) = HFloat(266.2987484570526), y2_4(x) = HFloat(-13.486001518995927), y2_5(x) = HFloat(-0.02999714213517951), y2_6(x) = HFloat(-22.43800529304631), y2_7(x) = HFloat(13.486001518918538), y2_8(x) = HFloat(-94.16612630071822), y2_9(x) = HFloat(785.015722665203), y2_10(x) = HFloat(-6775.009356439386), y2_11(x) = HFloat(-11.206236970774473), y2_12(x) = HFloat(77.76745611873598), y2_13(x) = HFloat(-642.9939052968533), y2_14(x) = HFloat(4.689768932388945e-16), y2_15(x) = HFloat(-0.0187710269573956), y2_16(x) = HFloat(0.0020861962139203305), y2_17(x) = HFloat(-5.312274942589151e-5), y2_18(x) = HFloat(-0.005830316541460037), y2_19(x) = HFloat(-1.227143374954293e-4), y2_20(x) = HFloat(-0.0011198874012043685), y2_21(x) = HFloat(7.602496797322661e-5), y2_22(x) = HFloat(-2.932591600936832e-6), y2_23(x) = HFloat(5.109018589948548e-7), y2_24(x) = HFloat(-0.0066417834273388736), y2_25(x) = HFloat(1.1855664085790522e-4), y2_26(x) = HFloat(-6.471131431633229e-6), y3_1(x) = HFloat(0.0), y3_2(x) = HFloat(0.0), y3_3(x) = HFloat(0.0), y3_4(x) = HFloat(0.0), y3_5(x) = HFloat(0.0), y3_6(x) = HFloat(0.0), y3_7(x) = HFloat(0.0), y3_8(x) = HFloat(0.0), y3_9(x) = HFloat(0.0), y3_10(x) = HFloat(0.0), y3_11(x) = HFloat(0.0), y3_12(x) = HFloat(0.0), y3_13(x) = HFloat(0.0), y3_14(x) = HFloat(-7.007769576850207e-6), y3_15(x) = HFloat(-0.004786448735922665), y3_16(x) = HFloat(4.689280929440655e-6), y3_17(x) = HFloat(-1.115140128648413e-4), y3_18(x) = HFloat(-0.004788926306328681), y3_19(x) = HFloat(4.716347003469115e-6), y3_20(x) = HFloat(-5.966958988267552e-5), y3_21(x) = HFloat(-4.710193875348028e-6), y3_22(x) = HFloat(3.723242564645378e-9), y3_23(x) = HFloat(1.7383541695404693e-11), y3_24(x) = HFloat(-0.004787700776648968), y3_25(x) = HFloat(-8.409570657167222e-8), y3_26(x) = HFloat(3.920653061680723e-13)]

 (1)

[x = 1., y1_1(x) = HFloat(0.0), y1_2(x) = HFloat(-3.420984483671309e-31), y1_3(x) = HFloat(-9.266873428961157e-30), y1_4(x) = HFloat(-2.8265796778082293), y1_5(x) = HFloat(4.585656258526918), y1_6(x) = HFloat(29.73356557713309), y1_7(x) = HFloat(2.8265796778082293), y1_8(x) = HFloat(1.6414688476893966), y1_9(x) = HFloat(78.13532261272604), y1_10(x) = HFloat(134.29458529790443), y1_11(x) = HFloat(-4.585656258526918), y1_12(x) = HFloat(-12.84587679453849), y1_13(x) = HFloat(-141.040214445319), y1_14(x) = HFloat(0.0), y1_15(x) = HFloat(-0.09999999999999995), y1_16(x) = HFloat(0.004086375487434791), y1_17(x) = HFloat(0.0), y1_18(x) = HFloat(0.0), y1_19(x) = HFloat(-2.458791784171863e-4), y1_20(x) = HFloat(-0.001957542333257033), y1_21(x) = HFloat(1.9956662771272542e-4), y1_22(x) = HFloat(-1.5079719874448682e-6), y1_23(x) = HFloat(3.3273624675095036e-8), y1_24(x) = HFloat(-7.164450668824305e-4), y1_25(x) = HFloat(1.1943543119401487e-5), y1_26(x) = HFloat(1.0902056689470582e-6), y1_27(x) = HFloat(-11.236234112909697), y2_1(x) = HFloat(7.740295426479526e-11), y2_2(x) = HFloat(11.236234112909697), y2_3(x) = HFloat(266.2987484570526), y2_4(x) = HFloat(-13.486001518995927), y2_5(x) = HFloat(-0.02999714213517951), y2_6(x) = HFloat(-22.43800529304631), y2_7(x) = HFloat(13.486001518918538), y2_8(x) = HFloat(-94.16612630071822), y2_9(x) = HFloat(785.015722665203), y2_10(x) = HFloat(-6775.009356439386), y2_11(x) = HFloat(-11.206236970774473), y2_12(x) = HFloat(77.76745611873598), y2_13(x) = HFloat(-642.9939052968533), y2_14(x) = HFloat(4.689768932388945e-16), y2_15(x) = HFloat(-0.0187710269573956), y2_16(x) = HFloat(0.0020861962139203305), y2_17(x) = HFloat(-5.312274942589151e-5), y2_18(x) = HFloat(-0.005830316541460037), y2_19(x) = HFloat(-1.227143374954293e-4), y2_20(x) = HFloat(-0.0011198874012043685), y2_21(x) = HFloat(7.602496797322661e-5), y2_22(x) = HFloat(-2.932591600936832e-6), y2_23(x) = HFloat(5.109018589948548e-7), y2_24(x) = HFloat(-0.0066417834273388736), y2_25(x) = HFloat(1.1855664085790522e-4), y2_26(x) = HFloat(-6.471131431633229e-6), y3_1(x) = HFloat(0.0), y3_2(x) = HFloat(0.0), y3_3(x) = HFloat(0.0), y3_4(x) = HFloat(0.0), y3_5(x) = HFloat(0.0), y3_6(x) = HFloat(0.0), y3_7(x) = HFloat(0.0), y3_8(x) = HFloat(0.0), y3_9(x) = HFloat(0.0), y3_10(x) = HFloat(0.0), y3_11(x) = HFloat(0.0), y3_12(x) = HFloat(0.0), y3_13(x) = HFloat(0.0), y3_14(x) = HFloat(-7.007769576850207e-6), y3_15(x) = HFloat(-0.004786448735922665), y3_16(x) = HFloat(4.689280929440655e-6), y3_17(x) = HFloat(-1.115140128648413e-4), y3_18(x) = HFloat(-0.004788926306328681), y3_19(x) = HFloat(4.716347003469115e-6), y3_20(x) = HFloat(-5.966958988267552e-5), y3_21(x) = HFloat(-4.710193875348028e-6), y3_22(x) = HFloat(3.723242564645378e-9), y3_23(x) = HFloat(1.7383541695404693e-11), y3_24(x) = HFloat(-0.004787700776648968), y3_25(x) = HFloat(-8.409570657167222e-8), y3_26(x) = HFloat(3.920653061680723e-13)]
``

 (2)

solnumeric2 := dsolve(eval(odesys, lambda = .2), numeric, approxsoln = solnumeric1)

solnumeric2(1)

[x = 1., y1_1(x) = HFloat(0.0), y1_2(x) = HFloat(-3.030574676828747e-31), y1_3(x) = HFloat(-8.154092896420496e-30), y1_4(x) = HFloat(-2.8265796778082306), y1_5(x) = HFloat(4.585656258526919), y1_6(x) = HFloat(29.733565577133103), y1_7(x) = HFloat(2.8265796778082306), y1_8(x) = HFloat(1.6414688476894284), y1_9(x) = HFloat(78.13532261272591), y1_10(x) = HFloat(134.2945852979062), y1_11(x) = HFloat(-4.585656258526919), y1_12(x) = HFloat(-12.845876794538517), y1_13(x) = HFloat(-141.04021444531904), y1_14(x) = HFloat(0.0), y1_15(x) = HFloat(-0.09999999999999996), y1_16(x) = HFloat(0.004086375487434796), y1_17(x) = HFloat(0.0), y1_18(x) = HFloat(0.0), y1_19(x) = HFloat(-2.458791784171864e-4), y1_20(x) = HFloat(-0.0019575423332570297), y1_21(x) = HFloat(1.9956662771272528e-4), y1_22(x) = HFloat(-1.50797198744487e-6), y1_23(x) = HFloat(3.3273624675095235e-8), y1_24(x) = HFloat(-7.164450668824302e-4), y1_25(x) = HFloat(1.1943543119401463e-5), y1_26(x) = HFloat(1.0902056689470559e-6), y1_27(x) = HFloat(-11.2362341129097), y2_1(x) = HFloat(7.741244519455937e-11), y2_2(x) = HFloat(11.2362341129097), y2_3(x) = HFloat(266.2987484570527), y2_4(x) = HFloat(-13.486001518995932), y2_5(x) = HFloat(-0.029997142135179216), y2_6(x) = HFloat(-22.438005293046288), y2_7(x) = HFloat(13.486001518918526), y2_8(x) = HFloat(-94.16612630071823), y2_9(x) = HFloat(785.0157226652033), y2_10(x) = HFloat(-6775.0093564393865), y2_11(x) = HFloat(-11.206236970774471), y2_12(x) = HFloat(77.76745611873608), y2_13(x) = HFloat(-642.9939052968541), y2_14(x) = HFloat(4.690343978496297e-16), y2_15(x) = HFloat(-0.018771026957395617), y2_16(x) = HFloat(0.002086196213920331), y2_17(x) = HFloat(-5.312274942589153e-5), y2_18(x) = HFloat(-0.005830316541460037), y2_19(x) = HFloat(-1.2271433749542907e-4), y2_20(x) = HFloat(-0.001119887401204369), y2_21(x) = HFloat(7.602496797322656e-5), y2_22(x) = HFloat(-2.9325916009368313e-6), y2_23(x) = HFloat(5.109018589948565e-7), y2_24(x) = HFloat(-0.006641783427338871), y2_25(x) = HFloat(1.1855664085790523e-4), y2_26(x) = HFloat(-6.471131431633234e-6), y3_1(x) = HFloat(0.0), y3_2(x) = HFloat(0.0), y3_3(x) = HFloat(0.0), y3_4(x) = HFloat(0.0), y3_5(x) = HFloat(0.0), y3_6(x) = HFloat(0.0), y3_7(x) = HFloat(0.0), y3_8(x) = HFloat(0.0), y3_9(x) = HFloat(0.0), y3_10(x) = HFloat(0.0), y3_11(x) = HFloat(0.0), y3_12(x) = HFloat(0.0), y3_13(x) = HFloat(0.0), y3_14(x) = HFloat(-7.007769576850145e-6), y3_15(x) = HFloat(-0.00478644873592267), y3_16(x) = HFloat(4.6892809294406756e-6), y3_17(x) = HFloat(-1.115140128648413e-4), y3_18(x) = HFloat(-0.004788926306328684), y3_19(x) = HFloat(4.716347003469143e-6), y3_20(x) = HFloat(-5.966958988267562e-5), y3_21(x) = HFloat(-4.710193875348042e-6), y3_22(x) = HFloat(3.7232425646457026e-9), y3_23(x) = HFloat(1.7383541695557647e-11), y3_24(x) = HFloat(-0.004787700776648964), y3_25(x) = HFloat(-8.409570657166629e-8), y3_26(x) = HFloat(3.9206530681961753e-13)]

 (3)

NULL

``


 

Download toughODE.mw

 

How many ways do you want?...

tomleslie 13876
September 24 2020
1 1

Four possible alternatives are shown in the attached

restart

with(Student:-VectorCalculus)

BasisFormat(true)

v := VectorField(`<,>`(x, y, z))

Vector(3, {(1) = x, (2) = y, (3) = z})

 (1)

ScalarPotential(v)

(1/2)*x^2+(1/2)*y^2+(1/2)*z^2

 (2)

 

Thats easy to use the ScalarPotential for a manual exercise

Example 2: finding a Potential function

Show that F = [exp(x)*cos(y)+y*z, x*z-exp(x)*sin(y), x*y+z]

F = [exp(x)*cos(y)+y*z, x*z-exp(x)*sin(y), x*y+z]

 (3)

 is conversative and find a potential function for it.

F := VectorField(`<,>`(exp(x)*cos(y)+y*z, x*z-exp(x)*sin(y), x*y+z))

Vector(3, {(1) = exp(x)*cos(y)+y*z, (2) = x*z-exp(x)*sin(y), (3) = x*y+z})

 (4)

"(->)"

y*z*x+exp(x)*cos(y)+(1/2)*z^2

 (5)

ScalarPotential(F); pdsolve(F = Gradient(SP(x, y, z))); pdsolve(F = Nabla(SP(x, y, z))); eqs := [seq(F[i] = Nabla(SP(x, y, z))[i], i = 1 .. 3)]; pdsolve(eqs)

{SP(x, y, z) = (1/2)*z^2+y*z*x+exp(x)*cos(y)+_C1}

 (6)

NULL

``


 

Download SP2.mw

Clarify...

tomleslie 13876
September 23 2020
1 3

If I open the help page DeepLearning:-Classify as a worksheet and then execcute it repeatedly - it always *seems* to work.

See the attached (which for some reason won't display inline here,  not sure why??)

So I think you will have to provide the worksheet which exhibits the failure (Use the big green  up-arrow in the Mapleprimes toolbar)

BTW I'm using 64-bit Win7

Download DLClass.mw

Basically...

tomleslie 13876
September 23 2020
1 1

You just need to perform computations to much higher accuracy - in the attached I used Digits=50 (whihc may(?) be overkill

It is also probably a good idea to use RootFinding:-Nextzero() rather than RootFinding:-Analytic(), since you are only searching for roots along the real line.

It is also much simpler to use the LinearAlgebra:-GenerateMatrix() command to construct the  matrix whose determinant is needed to obtain the roots.

See the attached

> 

restart; Digits := 50; interface(displayprecision = 10); with(LinearAlgebra); with(RootFinding)
> 

Phi := unapply(C1*cos(beta*x)+C2*sin(beta*x)+C3*cosh(beta*x)+C4*sinh(beta*x), x); bc1 := EI*((D@@2)(Phi))(0)-Kr*(D(Phi))(0); bc2 := EI*((D@@3)(Phi))(0)+Kt*Phi(0); bc3 := ((D@@2)(Phi))(L)/beta^2; bc4 := ((D@@3)(Phi))(L)/beta^3; M := GenerateMatrix([bc1, bc2, bc3, bc4], [C1, C2, C3, C4])[1]; detM := eval(simplify(Determinant(M)), [EI = 20*10^9*7243.32, L = 121, Kt = 4.5*10^9, Kr = 5.7*10^9, rho = 400, A = 30]); rts := Array(0 .. 5); resids := Array(1 .. 5); for i to 5 do rts[i] := NextZero(unapply(detM, beta), rts[i-1]); resids[i] := eval(detM, beta = rts[i]) end do; rts[1 .. 5]; resids[1 .. 5]

bc1 := EI*(-C1*beta^2+C3*beta^2)-Kr*(C2*beta+C4*beta)

  

bc2 := EI*(-C2*beta^3+C4*beta^3)+Kt*(C1+C3)

  

bc3 := (-C1*beta^2*cos(beta*L)-C2*beta^2*sin(beta*L)+C3*beta^2*cosh(beta*L)+C4*beta^2*sinh(beta*L))/beta^2

  

bc4 := (C1*beta^3*sin(beta*L)-C2*beta^3*cos(beta*L)+C3*beta^3*sinh(beta*L)+C4*beta^3*cosh(beta*L))/beta^3

  

Matrix(4, 4, {(1, 1) = -EI*beta^2, (1, 2) = -Kr*beta, (1, 3) = EI*beta^2, (1, 4) = -Kr*beta, (2, 1) = Kt, (2, 2) = -EI*beta^3, (2, 3) = Kt, (2, 4) = EI*beta^3, (3, 1) = -cos(beta*L), (3, 2) = -sin(beta*L), (3, 3) = cosh(beta*L), (3, 4) = sinh(beta*L), (4, 1) = sin(beta*L), (4, 2) = -cos(beta*L), (4, 3) = sinh(beta*L), (4, 4) = cosh(beta*L)})

  

Vector[row](5, {(1) = 0.28561117575040674433472658728587202863832735961218e-2, (2) = 0.27953120389950134576687262910544157423896370383064e-1, (3) = 0.43550314889288492708146739421724924920212998363216e-1, (4) = 0.65908358098826847376941652358667073390929887695958e-1, (5) = 0.91228672905353004995231758761707573469828658718856e-1})

  

Array(%id = 18446744074366223774)

 (1)
> 

``
> 

 


 

Download roots.mw

Very doubtful boundary condition...

tomleslie 13876
September 22 2020
1 0

you might want to consider the attached very carefully

> 

  restart:
#
# Define some parameters
#
   params:=[phi1 = .1, phi2 = .1, rhos1 = 2720, rhos2 = 2810,
            rhosf = 997.1, khnf = 1.083061737, kf = .613,
            cp1 = 893, cp2 = 960, cpf = 4179, Pr = 6.2,
            knf = .8154646474, S=0.5, R=0.5, gamma=0.5
          ]:
#
# Define the ODES
#
  ODES := (diff(f(eta), `$`(eta, 4)))/((1-phi1)^2.5*(1-phi2)^2.5*((1-phi2)*(1-phi1+phi1*rhos1/rhosf)+phi2*rhos2/rhosf))+S*(f(eta)*(diff(f(eta), `$`(eta, 3)))-3*(diff(f(eta), `$`(eta, 2)))-eta*(diff(f(eta), `$`(eta, 3)))-(diff(f(eta), eta))*(diff(f(eta), `$`(eta, 2)))) = 0,
          (khnf/kf+(4/3)*R)*(diff(theta(eta), `$`(eta, 2)))/((1-phi2)*(1-phi1+phi1*rhos1*cp1/(rhosf*cpf))+phi2*rhos2*cp2/(rhosf*cpf))+S*Pr*(f(eta)*(diff(theta(eta), eta))-eta*(diff(theta(eta), eta))-gamma*(eta^2*(diff(theta(eta), `$`(eta, 2)))-2*eta*f(eta)*(diff(theta(eta), `$`(eta, 2)))-eta*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)*(diff(f(eta), eta))*(diff(theta(eta), eta))+f(eta)^2*(diff(theta(eta), `$`(eta, 2))))) = 0:
#
# Define the boundary conditions. Notice that in the OPs
# original definition ((D^2)(f))(0) = 0 is a a bit dumb.
# Apart from anything else it could be more simply written
# written as D(f)(0) = 0
#
# Can only suggest that OP reads the output of these three
# commands very carefully and decide which is required

  BCs0:= f(0) = 0, ((D^2)(f))(0) = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1;
  BCs1:= f(0) = 0, D(f)(0) = 0, D(theta)(0) = 0, f(1) = 0, D(f)(1) = 0, theta(1) = 1;
  BCs2:= f(0) = 0, D[1,1](f)(0) = 0, D(theta)(0) = 0, f(1) = 0, D(f)(1) = 0, theta(1) = 1;

f(0) = 0, (D(f))(0)^2 = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1

  

f(0) = 0, (D(f))(0) = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1

  

f(0) = 0, ((D@@2)(f))(0) = 0, (D(theta))(0) = 0, f(1) = 0, (D(f))(1) = 0, theta(1) = 1

 (1)
> 

#
# Since I don't think the first set of BCs above is plausible
# because of the existence of ((D^2)(f))(0) = 0, just solve the
# system with the second and third set of boundary conditions.
# Neither produces very interesting solutions
#
  sol1:=dsolve( [eval([ODES], params)[], BCs1], numeric);
  sol2:=dsolve( [eval([ODES], params)[], BCs2], numeric);
  plots:-odeplot(sol1, [[eta, f(eta)], [eta,theta(eta)]], eta=0..1, axes=boxed, color=[red,blue]);
  plots:-odeplot(sol2, [[eta, f(eta)], [eta,theta(eta)]], eta=0..1, axes=boxed, color=[red,blue]);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = 1.0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = 1.0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = 1.0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = 1.0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = 1.0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = 1.0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = 1.0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = 1.0, (8, 6) = .0}, datatype = float[8], order = C_order); YP := Matrix(8, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [6, 8, [f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(6, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(6, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = yout[i], i = 1 .. 6)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [6, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(6, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0.}); `dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 6)] end proc, (2) = Array(0..0, {}), (3) = [eta, f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [eta = res[1], seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = res[i+1], i = 1 .. 6)] catch: error end try end proc

  

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = 1.0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = 1.0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = 1.0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = 1.0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = 1.0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = 1.0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = 1.0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = 1.0, (8, 6) = .0}, datatype = float[8], order = C_order); YP := Matrix(8, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = .14285714285714277, (3) = .28571428571428553, (4) = .4285714285714284, (5) = .5714285714285713, (6) = .7142857142857142, (7) = .857142857142857, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 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, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [6, 8, [f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(6, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(6, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 6, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = yout[i], i = 1 .. 6)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return HFloat(-0.0) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [6, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(6, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(6, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0.}); `dsolve/numeric/hermite`(8, 6, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 6)] end proc, (2) = Array(0..0, {}), (3) = [eta, f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [eta = res[1], seq('[f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta), diff(diff(diff(f(eta), eta), eta), eta), theta(eta), diff(theta(eta), eta)]'[i] = res[i+1], i = 1 .. 6)] catch: error end try end proc

  

 

 
> 

 

> 

 

Download odeProb.mw

confusion...

tomleslie 13876
September 20 2020
0 6

over 2D-input and 2D-math formats

Is the attached what you expected (For some reason this wont display inline  maybe the animations make the file size too big?)

anim.mw

Something like...

tomleslie 13876
September 20 2020
0 1

@AHSAN 

the attached, maybe?

> 

restart
> 

A := beta*(-60*lambda*(1+(1/2)*x^2)^4+(-63*k^3+99*k^2-99*k+63)*(1+(1/2)*x^2)^3-288*Q*(k^2-(11/8)*k+1)*(1+(1/2)*x^2)^2-486*Q^2*(k-1)*(1+(1/2)*x^2)-324*Q^3)/(20*(1+(1/2)*x^2)^7)-(3*(k*(1+(1/2)*x^2)+2*Q-1-(1/2)*x^2))/(2*(1+(1/2)*x^2)^3)
> 

opts := [[beta = .1, Q = .5516, lambda = -0.396e-1, k = .1], [beta = .2, Q = .4903, lambda = -0.411e-1, k = .2], [beta = .3, Q = .4290, lambda = -0.419e-1, k = .3], [beta = .4, Q = .3677, lambda = -0.417e-1, k = .4]]
> 

yVals := [-.1, -.12, -.14, -.16]; N := plot([seq(eval(A, opts[j]), j = 1 .. 4)], x = -2 .. 2, color = [black, red, blue, green], thickness = 2); B := plot(yVals, x = -1 .. -.7, color = [black, red, blue, green], thickness = 3); T := plots:-textplot([seq([-1.1, yVals[j], typeset(k = rhs(opts[j, 3]), ", ", Q = rhs(opts[j, 2]))], j = 1 .. 4)], align = left)
> 

plots:-display(N, B, T, size = [900, 700], gridlines = false)

 
> 

 


 

Download legend2.mw

Maybe you woud prefer the output from if...

tomleslie 13876
September 16 2020
0 3

@user-JG 

as in

> 

ifactor(725);
ifactor(1125);
ifactor(2048);

``(5)^2*``(29)

  

``(3)^2*``(5)^3

  

``(2)^11

 (1)
> 

 


 

Download fac.mw

First 65 66 67 68 69 70 71 Last Page 67 of 207