Advanced Search

dharr

Dr. David Harrington

8195 Reputation

22 Badges

20 years, 334 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

Double sum...

dharr 8195
March 24 2025
2 10

The final answer is

sum(sum(1/(m^2 - n^2), m = 1 .. n - 1)
 + sum(1/(m^2 - n^2), m = n + 1 .. infinity), n = 1 .. infinity)

giving Pi^2/8. Perhaps obvious, but I got it by hacking around as below, which might or might not be instructive.

Summation.mw

comments...

dharr 8195
March 23 2025
1 2

A large section of 2D code can be harder to delete or copy (or paste as 2D), since internally it has a longer format. If you want things to be faster to process, use ctrl-M at the prompt before pasting it in from another worksheet, so it is converted to 1D for later faster use.

The first plot is just the textplot T. To not show it, end the T:= statement with a colon (":").

Check denominators...

dharr 8195
March 21 2025
2 0

For example you have a denominator 

-3*(k[1] + k[2])^2*alpha + beta*(l[1] - l[2])^2

but k[1]=-k[2] and l[1]=l[2] so this evaluates to zero.

via generating functions...

dharr 8195
March 21 2025
1 0

Since Maple has difficulty with symbolic k, one can work it out for many k and then find a potential answer using generating functions.

Download generating_functions.mw

without the cat...

dharr 8195
March 19 2025
1 0

Possibly a version issue, but for me in Maple 2024.2 it works without cat. Or perhaps you want something else?

GenerateSimilar_Test.mw

removeable singularity...

dharr 8195
March 12 2025
1 0

This solution can be found by using other equivalent initial conditions, presumably even ones close to but not exactly at t=0. But y(0)=1 is more sensible than y(0)=0.

> 

restart;
> 

ode:=(1/(1+t^2)-y(t)^2)-(2*t*y(t))*diff(y(t),t);
IC:=y(0)=0;

1/(t^2+1)-y(t)^2-2*t*y(t)*(diff(y(t), t))

y(0) = 0

> 

sol:=dsolve([ode,IC],numeric);

Error, (in dsolve/numeric/checksing) ode system has a removable singularity at t=0. Initial data is restricted to one of the following 2 options: {{y(t) = -1}, {y(t) = 1}}

Solution seems to be defined nicely everywhere but zero

> 

z:=sqrt(arctan(t)/t);
limit(z,t=0);
plot(z,t=-2..2);

(arctan(t)/t)^(1/2)

1

> 

IC2:=y(1)=eval(z,t=1);

y(1) = (1/4)*4^(1/2)*Pi^(1/2)

> 

solnum:=dsolve([ode,IC2],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 .. 28, [( 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..65, {(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, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = 1.0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = 1.0, (6) = .30916520090870503, (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) = .88622692545276}, 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) = .88622692545276}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = -.16101867095250255}, 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] = y(t)]`; YP[1] := (1/2)*(1/(X^2+1)-Y[1]^2)/(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(t)]`; YP[1] := (1/2)*(1/(X^2+1)-Y[1]^2)/(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..1, {(1) = 1.}); _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 elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) 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) = [t, y(t)], (4) = []}); _vars := _dat[3]; _pars := map(lhs, _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

> 

plots:-odeplot(solnum,[t,y(t)],0..2);

Warning, cannot evaluate the solution further left of .44217005e-6, probably a singularity

> 

sol2:=dsolve([ode,IC2]);

y(t) = (t*arctan(t))^(1/2)/t

Despite the message from dsolve/numeric, other initial condition locations also work.

> 

IC3:=y(1/2)=eval(z, t=1/2);
sol3:=dsolve([ode,IC3]);

y(1/2) = 2^(1/2)*arctan(1/2)^(1/2)

y(t) = (t*arctan(t))^(1/2)/t

> 

IC4:=y(0)=1;
sol4:=dsolve([ode,IC4]);

y(0) = 1

> 

NULL

Download dsolve.mw

simplify with side relations...

dharr 8195
March 11 2025
2 2
> 

restart
> 

with(VectorCalculus)
> 

with(DETools)
> 

with(DynamicSystems)
> 

_local(I)

I

Warning, The imaginary unit, I, has been renamed _I

> 

dS := -P*S*alpha-S*mu+Lambda; dE := alpha*S*P-(-T*eta+1)*beta*E-theta*E-mu*E; dI := (-T*eta+1)*beta*E-delta*I-gamma*I-mu*I; dR := E*theta+I*gamma-R*mu; dP := -P*T*sigma+I*xi-P*tau; dT := r*T*(1-T/K)-phi*T

-P*S*alpha-S*mu+Lambda

alpha*S*P-(-T*eta+1)*beta*E-theta*E-mu*E

(-T*eta+1)*beta*E-delta*I-gamma*I-mu*I

E*theta+I*gamma-R*mu

-P*T*sigma+I*xi-P*tau

r*T*(1-T/K)-phi*T

> 

SEkuil := solve({dE, dI, dP, dR, dS, dT}, {E, I, P, R, S, T})
> 

SEkuil_End := SEkuil[4]
> 

E1 := eval(E, SEkuil_End)

(K^2*beta*delta*eta*mu*phi^2*sigma-2*K^2*beta*delta*eta*mu*phi*r*sigma+K^2*beta*delta*eta*mu*r^2*sigma+K^2*beta*eta*gamma*mu*phi^2*sigma-2*K^2*beta*eta*gamma*mu*phi*r*sigma+K^2*beta*eta*gamma*mu*r^2*sigma+K^2*beta*eta*mu^2*phi^2*sigma-2*K^2*beta*eta*mu^2*phi*r*sigma+K^2*beta*eta*mu^2*r^2*sigma+K*Lambda*alpha*beta*eta*phi*r*xi-K*Lambda*alpha*beta*eta*r^2*xi-K*beta*delta*eta*mu*phi*r*tau+K*beta*delta*eta*mu*r^2*tau-K*beta*eta*gamma*mu*phi*r*tau+K*beta*eta*gamma*mu*r^2*tau-K*beta*eta*mu^2*phi*r*tau+K*beta*eta*mu^2*r^2*tau+K*beta*delta*mu*phi*r*sigma-K*beta*delta*mu*r^2*sigma+K*beta*gamma*mu*phi*r*sigma-K*beta*gamma*mu*r^2*sigma+K*beta*mu^2*phi*r*sigma-K*beta*mu^2*r^2*sigma+K*delta*mu^2*phi*r*sigma-K*delta*mu^2*r^2*sigma+K*delta*mu*phi*r*sigma*theta-K*delta*mu*r^2*sigma*theta+K*gamma*mu^2*phi*r*sigma-K*gamma*mu^2*r^2*sigma+K*gamma*mu*phi*r*sigma*theta-K*gamma*mu*r^2*sigma*theta+K*mu^3*phi*r*sigma-K*mu^3*r^2*sigma+K*mu^2*phi*r*sigma*theta-K*mu^2*r^2*sigma*theta+Lambda*alpha*beta*r^2*xi-beta*delta*mu*r^2*tau-beta*gamma*mu*r^2*tau-beta*mu^2*r^2*tau-delta*mu^2*r^2*tau-delta*mu*r^2*tau*theta-gamma*mu^2*r^2*tau-gamma*mu*r^2*tau*theta-mu^3*r^2*tau-mu^2*r^2*tau*theta)/((K*eta*phi-K*eta*r+r)*xi*beta*(K*beta*eta*phi-K*beta*eta*r+beta*r+mu*r+r*theta)*alpha)

We want the above solution in terms of the following expression.

> 

R0 := Lambda*alpha*beta*r*xi*(-K*eta*phi+K*eta*r-r)/(mu*(K*phi*sigma-K*r*sigma-r*tau)*(K*beta*eta*phi-K*beta*eta*r+beta*r+mu*r+r*theta)*(delta+gamma+mu))

Lambda*alpha*beta*r*xi*(-K*eta*phi+K*eta*r-r)/(mu*(K*phi*sigma-K*r*sigma-r*tau)*(K*beta*eta*phi-K*beta*eta*r+beta*r+mu*r+r*theta)*(delta+gamma+mu))

We can use simplify with side relations but we want the above relationship in polynomial form. This can be achieved by multiplying each side by the denominator, which amounts to requiring the following to be zero

> 

siderel := numer(normal(R__0-R0))

K^2*R__0*beta*delta*eta*mu*phi^2*sigma-2*K^2*R__0*beta*delta*eta*mu*phi*r*sigma+K^2*R__0*beta*delta*eta*mu*r^2*sigma+K^2*R__0*beta*eta*gamma*mu*phi^2*sigma-2*K^2*R__0*beta*eta*gamma*mu*phi*r*sigma+K^2*R__0*beta*eta*gamma*mu*r^2*sigma+K^2*R__0*beta*eta*mu^2*phi^2*sigma-2*K^2*R__0*beta*eta*mu^2*phi*r*sigma+K^2*R__0*beta*eta*mu^2*r^2*sigma-K*R__0*beta*delta*eta*mu*phi*r*tau+K*R__0*beta*delta*eta*mu*r^2*tau-K*R__0*beta*eta*gamma*mu*phi*r*tau+K*R__0*beta*eta*gamma*mu*r^2*tau-K*R__0*beta*eta*mu^2*phi*r*tau+K*R__0*beta*eta*mu^2*r^2*tau+K*Lambda*alpha*beta*eta*phi*r*xi-K*Lambda*alpha*beta*eta*r^2*xi+K*R__0*beta*delta*mu*phi*r*sigma-K*R__0*beta*delta*mu*r^2*sigma+K*R__0*beta*gamma*mu*phi*r*sigma-K*R__0*beta*gamma*mu*r^2*sigma+K*R__0*beta*mu^2*phi*r*sigma-K*R__0*beta*mu^2*r^2*sigma+K*R__0*delta*mu^2*phi*r*sigma-K*R__0*delta*mu^2*r^2*sigma+K*R__0*delta*mu*phi*r*sigma*theta-K*R__0*delta*mu*r^2*sigma*theta+K*R__0*gamma*mu^2*phi*r*sigma-K*R__0*gamma*mu^2*r^2*sigma+K*R__0*gamma*mu*phi*r*sigma*theta-K*R__0*gamma*mu*r^2*sigma*theta+K*R__0*mu^3*phi*r*sigma-K*R__0*mu^3*r^2*sigma+K*R__0*mu^2*phi*r*sigma*theta-K*R__0*mu^2*r^2*sigma*theta-R__0*beta*delta*mu*r^2*tau-R__0*beta*gamma*mu*r^2*tau-R__0*beta*mu^2*r^2*tau-R__0*delta*mu^2*r^2*tau-R__0*delta*mu*r^2*tau*theta-R__0*gamma*mu^2*r^2*tau-R__0*gamma*mu*r^2*tau*theta-R__0*mu^3*r^2*tau-R__0*mu^2*r^2*tau*theta+Lambda*alpha*beta*r^2*xi

> 

simplify(E1, {siderel})

-(sigma*(r-phi)*K+r*tau)*(delta+gamma+mu)*(R__0-1)*mu/(xi*alpha*beta*(eta*(r-phi)*K-r))

> 

NULL

Download end.mw

a possibility...

dharr 8195
March 08 2025
1 3

I don't think there is a simple way to do this, but here is one way to make it convincing. Verifying the actual and your proposed formulas solve the equation is easier and could be a first step.

> 

restart
> 

eqn1 := sin(x) = 1/2

sin(x) = 1/2

Solve and then replace the variable starting with _Z by k.

> 

ans1, ans2 := solve(eqn1, allsolutions); ans1 := eval(ans1, indets(ans1, suffixed(_Z))[] = k); ans2 := eval(ans2, indets(ans2, suffixed(_Z))[] = k)

(1/6)*Pi+2*Pi*_Z1, (5/6)*Pi+2*Pi*_Z1

(1/6)*Pi+2*Pi*k

(5/6)*Pi+2*Pi*k

> 

ans3 := k*Pi+(1/6)*(-1)^k*Pi

Pi*k+(1/6)*(-1)^k*Pi

Check each of them individually satisfies the equation (could also apply evalb or is to get true).

> 

`assuming`([simplify(eval(eqn1, x = ans1))], [k::integer]); `assuming`([simplify(eval(eqn1, x = ans2))], [k::integer]); `assuming`([simplify(eval(eqn1, x = ans3))], [k::integer])

1/2 = 1/2

1/2 = 1/2

1/2 = 1/2

Of course this doesn't show ans1 and ans2 cover the same values as ans3. We can generate some possibilities to see if this is true, but you have to try different ranges of k to see that they match

> 

seq(ans3, k = -5 .. 5); seq(ans1, k = -2 .. 2); seq(ans2, k = -3 .. 2)

-(31/6)*Pi, -(23/6)*Pi, -(19/6)*Pi, -(11/6)*Pi, -(7/6)*Pi, (1/6)*Pi, (5/6)*Pi, (13/6)*Pi, (17/6)*Pi, (25/6)*Pi, (29/6)*Pi

-(23/6)*Pi, -(11/6)*Pi, (1/6)*Pi, (13/6)*Pi, (25/6)*Pi

-(31/6)*Pi, -(19/6)*Pi, -(7/6)*Pi, (5/6)*Pi, (17/6)*Pi, (29/6)*Pi

So this works

> 

{seq(ans1, k = -2 .. 2), seq(ans2, k = -3 .. 2)} = {seq(ans3, k = -5 .. 5)}; evalb(%)

{-(31/6)*Pi, -(23/6)*Pi, -(19/6)*Pi, -(11/6)*Pi, -(7/6)*Pi, (1/6)*Pi, (5/6)*Pi, (13/6)*Pi, (17/6)*Pi, (25/6)*Pi, (29/6)*Pi} = {-(31/6)*Pi, -(23/6)*Pi, -(19/6)*Pi, -(11/6)*Pi, -(7/6)*Pi, (1/6)*Pi, (5/6)*Pi, (13/6)*Pi, (17/6)*Pi, (25/6)*Pi, (29/6)*Pi}

true

> 

NULL

Download sets.mw

convert(...,trigh)...

dharr 8195
March 07 2025
0 0

The answer to your second question is convert(Bans, trigh). But this can leave complex expressions as arguments, so you might want to also try convert(..., trig).

via products...

dharr 8195
March 06 2025
4 0

I more-or-less did it by hand here; I'm sure it could be streamlined.

> 

restart;
> 

s1 := arctan(2/n^2);
s2 := expand(convert(s1, ln));

arctan(2/n^2)

((1/2)*I)*ln(1-(2*I)/n^2)-((1/2)*I)*ln(1+(2*I)/n^2)

Take exp so we can do products instead of sums; find the pieces

> 

p:=simplify(exp(s2),exp);
p1,e1:=op(op(1,p));
p2,e2:=op(op(2,p));

(1-(2*I)/n^2)^((1/2)*I)*(1+(2*I)/n^2)^(-(1/2)*I)

1-(2*I)/n^2, (1/2)*I

1+(2*I)/n^2, -(1/2)*I

We can do the products

> 

P1:=product(p1, n=1.. infinity);
P2:=product(p2, n=1.. infinity);

(-1/2+(1/2)*I)*sin((2+I)*Pi)/Pi

(-1/2-(1/2)*I)*sin((2-I)*Pi)/Pi

> 

Take the ln to reverse the expNULL
> 

S:=simplify(e1*ln(P1)+e2*ln(P2));

(3/4)*Pi

> 

 

Download sum.mw

Some tips...

dharr 8195
March 03 2025
1 0

Here's how I would do it, keeping close to how you did it. When using solve, if you get a warning about assumptions not used, add the useassumptions option.

restart

Levine - Ch. 2.4 - Particle in a Rectangular Well

 

In this problem, E < V__0 and m > 0.

assume(E < V__0, m > 0, `&hbar;` > 0)

NULL

Define

s__1 := sqrt(2*m*(V__0-E))/`&hbar;` = 2^(1/2)*(m*(V__0-E))^(1/2)/`&hbar;`NULL

s__2 := -s__1 = -2^(1/2)*(m*(V__0-E))^(1/2)/`&hbar;`NULL

s := sqrt(2*m*E)/`&hbar;`

2^(1/2)*(m*E)^(1/2)/`&hbar;`

`&psi;__1`, `&psi;__2` and `&psi;__3` are wavefunctions.

`&psi;__1` := proc (x) options operator, arrow; C*exp('s__1'*x) end proc = proc (x) options operator, arrow; C*exp('s__1'*x) end procNULL

`&psi;__2` := proc (x) options operator, arrow; A*cos('s'*x)+B*sin('s'*x) end proc = proc (x) options operator, arrow; A*cos('s'*x)+B*sin('s'*x) end procNULL

`&psi;__3` := proc (x) options operator, arrow; G*exp('s__2'*x) end proc = proc (x) options operator, arrow; G*exp('s__2'*x) end procNULL

``

As can be seen above, there are four unknown constants, "A,B,C,"and G.

 

We can obtain values for these unknowns by imposing boundary conditions.

 

Some of the boundary conditions involve the first derivatives of the wavefunctions.

`&psi;__1,d` := D(`&psi;__1`) = proc (z) options operator, arrow; C*s__1*exp(s__1*z) end procNULL

`&psi;__2,d` := D(`&psi;__2`) = proc (z) options operator, arrow; -A*s*sin(s*z)+B*s*cos(s*z) end procNULL

`&psi;__3,d` := D(`&psi;__3`) = proc (z) options operator, arrow; G*s__2*exp(s__2*z) end procNULL

NULL

The boundary conditions are

NULL

`&psi;__1`(0) = `&psi;__2`(0)*`&psi;__2`(l) and `&psi;__2`(0)*`&psi;__2`(l) = `&psi;__3`(l)*(D(`&psi;__1`))(0) and `&psi;__3`(l)*(D(`&psi;__1`))(0) = (D(`&psi;__2`))(0)*(D(`&psi;__2`))(l) and (D(`&psi;__2`))(0)*(D(`&psi;__2`))(l) = (D(`&psi;__3`))(l)

NULL

My question is how to find the constants in Maple.

NULL``

Solving the first boundary condition is easy.

expr1 := solve(`&psi;__1`(0) = `&psi;__2`(0), {C}, useassumptions) = {C = A}NULL

NULL

Next, I solve the third boundary condition and sub in the result from solving the first boundary condition. We get B in terms of A.

expr3 := simplify(eval(solve(`&psi;__1,d`(0) = `&psi;__2,d`(0), {B}), expr1)) = {B = A*(V__0-E)^(1/2)/E^(1/2)}NULL

``

Already here it is not clear to me why Maple does not cancel the m's - use simplify.

 

Next, consider the second boundary condition

expr2 := solve(`&psi;__2`(l) = `&psi;__3`(l), {G}, useassumptions) = {G = (A*cos(s*l)+B*sin(s*l))/exp(s__2*l)}NULL

 

We can obtain G in terms of just A by subbing in previously found relationships, as follows

expr2 := simplify(subs(expr3, expr2))

{G = exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)*A*(cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2)+(V__0-E)^(1/2)*sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`))/E^(1/2)}

NULL

Finally, consider the fourth boundary condition.

solve(`&psi;__2,d`(l) = `&psi;__3,d`(l), {G}, useassumptions)

{G = (m*E)^(1/2)*(A*sin(2^(1/2)*(m*E)^(1/2)*l/`&hbar;`)-B*cos(2^(1/2)*(m*E)^(1/2)*l/`&hbar;`))/((m*(V__0-E))^(1/2)*exp(-2^(1/2)*(m*(V__0-E))^(1/2)*l/`&hbar;`))}

expr4 := simplify(subs(expr3, solve(`&psi;__2,d`(l) = `&psi;__3,d`(l), {G}, useassumptions)))

{G = (-(V__0-E)^(1/2)*cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)+sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2))*A*exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)/(V__0-E)^(1/2)}

result := simplify(subs(expr2, expr4))

{exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)*A*(cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2)+(V__0-E)^(1/2)*sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`))/E^(1/2) = (-(V__0-E)^(1/2)*cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)+sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2))*A*exp(2^(1/2)*m^(1/2)*(V__0-E)^(1/2)*l/`&hbar;`)/(V__0-E)^(1/2)}

Subtract rhs-lhs, and divide through or multiply by the factors to cancel

result2 := simplify((rhs-lhs)(result[])*sqrt(V__0-E)*sqrt(E)/(exp(sqrt(2)*sqrt(m)*sqrt(V__0-E)*l/`&hbar;`)*A))

(2*E-V__0)*sin(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)-2*cos(2^(1/2)*m^(1/2)*E^(1/2)*l/`&hbar;`)*E^(1/2)*(V__0-E)^(1/2)

NULL

For example, I would like to replace sqrt(2*mE)/`&hbar;` with s. But s already has a value, so you need to use a different symbol or unassign it

s := 's'

s

Since this is essentially replacing all the ms

s = sqrt(2*m*E)/`&hbar;`; isolate(%, m); simplify(eval(result2, %))

(2*E-V__0)*sin((s^2/E)^(1/2)*E^(1/2)*l)-2*cos((s^2/E)^(1/2)*E^(1/2)*l)*E^(1/2)*(V__0-E)^(1/2)

 

This is a first task (that I have asked about in a separate question).

 

Then, I would want Maple to do the cancellations for me, but I have already used "simplify".

 

Finally, I think I would like to collect terms involving cos and sin.

 

Here is what happens when I try to get Maple to solve the equations for me in one go

 

solve({`&psi;__1`(0) = `&psi;__2`(0), `&psi;__1,d`(0) = `&psi;__2,d`(0), `&psi;__2`(l) = `&psi;__3`(l), `&psi;__2,d`(l) = `&psi;__3,d`(l)}, {A, B, C, G}, useassumptions)

{A = 0, B = 0, C = 0, G = 0}

NULL

 

NULL

Download solve_constants.mw

Error fixed; LambertW...

dharr 8195
February 28 2025
2 1

You did not change u(x,y,t) to u(x,y,z,t) in the code that extracts the linear part. After fixing that, there is no solution.
(The answer to the question about removing LambertW is that you can't).

remove2.mw

using indets...

dharr 8195
February 26 2025
2 2

Basically, using indets gives you the (more complicated) subexpressions that you can directly substitute. I did it twice so there are no exps or square roots left, but you can go further if you wish.

restart

with(LinearAlgebra)

 

Additional conditions:

c1 := `&omega;__1`+`&omega;__2` = 1; c2 := f__p+f__d1+f__d2 = 1

f__p+f__d1+f__d2 = 1

All non-negative values (but f__p and T must be positive for the denominators to make sense).NULL

NULL

F := Matrix([[-(f__d1/f__p+1)*`&lambda;__1`, -`&omega;__1`*f__d2*`&lambda;__2`/f__p, `&omega;__1`, 0], [-(`&omega;__2`*f__d1/(`&omega;__1`*f__p)-1)*`&lambda;__1`, -(`&omega;__2`*f__d2/f__p+1)*`&lambda;__2`, `&omega;__2`, 0], [0, 0, 0, 1/T], [0, 0, 0, 0]])

Matrix(%id = 36893490587109511524)

 

 

 

Make the following substitutions right away, to reduce the number of variables by 2

NULLc3 := {isolate(c1, `&omega;__1`), isolate(c2, f__d1)}

{f__d1 = 1-f__p-f__d2, omega__1 = 1-omega__2}

F2 := simplify(eval(F, c3))

Matrix(%id = 36893490587109494780)

expFT := MatrixExponential(F2, T)NULL

NULL

inds := indets(expFT)

{T, f__d2, f__p, lambda__1, lambda__2, omega__2, 1/((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2), ((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2), exp((1/2)*(-f__d2*lambda__2*omega__2+f__d2*lambda__1-f__p*lambda__2-lambda__1+((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2))*T/f__p), exp(-(1/2)*(f__d2*lambda__2*omega__2-f__d2*lambda__1+f__p*lambda__2+((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2)+lambda__1)*T/f__p)}

musubs := `~`[`=`](select(has, inds, exp), {mu[1], mu[2]})

{exp((1/2)*(-f__d2*lambda__2*omega__2+f__d2*lambda__1-f__p*lambda__2-lambda__1+((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2))*T/f__p) = mu[1], exp(-(1/2)*(f__d2*lambda__2*omega__2-f__d2*lambda__1+f__p*lambda__2+((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2)+lambda__1)*T/f__p) = mu[2]}

expFTmu := subs(musubs, expFT); inds2 := indets(expFTmu)

{T, f__d2, f__p, lambda__1, lambda__2, omega__2, mu[1], mu[2], 1/((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2), ((-1+f__d2)^2*lambda__1^2-2*(f__d2^2*omega__2+(f__p-omega__2)*f__d2+f__p)*lambda__2*lambda__1+lambda__2^2*(f__d2*omega__2+f__p)^2)^(1/2)}

``

Check last two are reciprocals; since they are, convert them to 1/P and P

simplify(inds2[-2]-1/inds2[-1]); expFTmuP := simplify(subs(inds2[-1] = P, inds2[-2] = 1/P, expFTmu))

Matrix(%id = 36893490587266132316)

``

NULL

Download expFT_compact.mw

P.S. If the matrix were diagonalizable, I could probably do something more intelligent via the eigenvectors. Is the last row really intended to be zero?

convert(..,radical)...

dharr 8195
February 25 2025
2 0

The normal way to deal with a RootOf if you want an explicit solution would be to convert it to radical form. In your case, it is the root(s) of a quadratic, so you can use convert(E[u], radical) to get an explicit solution in terms of the quadratic formula. If you want both roots, you can use allvalues(E[u]). Or better yet, you probably got that result from solve, so you could add the option explicit to the solve command and not generate the RootOf in the first place. (remove_RootOf is just giving you an equivalent equation to solve, so you are not getting ahead by that method).

In future please upload your worksheet using the green up-arrow in the Mapleprimes editor, select your file, click upload, and then insert link or insert contents.

One lump answer...

dharr 8195
February 21 2025
1 51

Here the derivation for the one lump case (sections in green). The two-lump case can be done similarly, I assume.

> 

restart
> 

with(PDEtools)
> 

undeclare(prime, quiet); declare(u(x, y, t), quiet); declare(f(x, y, t), quiet)
> 

pde := diff(u(x, y, t), t)-(diff(diff(u(x, y, t), `$`(x, 4))+5*u(x, y, t)*(diff(u(x, y, t), `$`(x, 2)))+(5/3)*u(x, y, t)^3+5*(diff(u(x, y, t), x, y)), x))-5*u(x, y, t)*(diff(u(x, y, t), y))+5*(int(diff(u(x, y, t), `$`(y, 2)), x))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

diff(u(x, y, t), t)-(diff(diff(diff(diff(diff(u(x, y, t), x), x), x), x), x))-5*(diff(u(x, y, t), x))*(diff(diff(u(x, y, t), x), x))-5*u(x, y, t)*(diff(diff(diff(u(x, y, t), x), x), x))-5*u(x, y, t)^2*(diff(u(x, y, t), x))-5*(diff(diff(diff(u(x, y, t), x), x), y))-5*u(x, y, t)*(diff(u(x, y, t), y))+5*(int(diff(diff(u(x, y, t), y), y), x))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

> 

pde_linear, pde_nonlinear := selectremove(proc (term) options operator, arrow; not has((eval(term, u(x, y, t) = a*u(x, y, t)))/a, a) end proc, expand(pde))

diff(u(x, y, t), t)-(diff(diff(diff(diff(diff(u(x, y, t), x), x), x), x), x))-5*(diff(diff(diff(u(x, y, t), x), x), y))+5*(int(diff(diff(u(x, y, t), y), y), x)), -5*(diff(u(x, y, t), x))*(diff(diff(u(x, y, t), x), x))-5*u(x, y, t)*(diff(diff(diff(u(x, y, t), x), x), x))-5*u(x, y, t)^2*(diff(u(x, y, t), x))-5*u(x, y, t)*(diff(u(x, y, t), y))-5*(diff(u(x, y, t), x))*(int(diff(u(x, y, t), y), x))

> 

thetai := t*w[i]+y*p[i]+x

t*w[i]+y*p[i]+x

> 

eqw := w[i] = -5*p[i]^2

w[i] = -5*p[i]^2

> 

Bij := proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

proc (i, j) options operator, arrow; (-6*p[i]-6*p[j])/(p[i]-p[j])^2 end proc

> 

NULL
> 

theta1 := normal(eval(eval(thetai, eqw), i = 1)); theta2 := normal(eval(eval(thetai, eqw), i = 2))

-5*t*p[1]^2+y*p[1]+x

-5*t*p[2]^2+y*p[2]+x

> 

eqf := f(x, y, t) = theta1*theta2+Bij(1, 2)

f(x, y, t) = (-5*t*p[1]^2+y*p[1]+x)*(-5*t*p[2]^2+y*p[2]+x)+(-6*p[1]-6*p[2])/(p[1]-p[2])^2

> 

eqfcomplex := collect(evalc(eval(eval(eqf, p[2] = conjugate(p[1])), p[1] = a+I*b)), t)

f(x, y, t) = ((-5*a^2+5*b^2)^2+100*a^2*b^2)*t^2+(2*(a*y+x)*(-5*a^2+5*b^2)-20*y*b^2*a)*t+(a*y+x)^2+y^2*b^2+3*a/b^2

Eqs 16,17

> 

eqp := {x = xp-(5*(a^2+b^2))*t, y = 10*a*t+yp}; eqfc2 := f(x, y, t) = simplify(eval(rhs(eqfcomplex), eqp))

f(x, y, t) = (b^4*yp^2+(a*yp+xp)^2*b^2+3*a)/b^2

> 

eq17 := u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2

u(x, y, t) = 6*(diff(diff(f(x, y, t), x), x))/f(x, y, t)-6*(diff(f(x, y, t), x))^2/f(x, y, t)^2

Eq 18

> 

eqt := u(x, y, t) = simplify(eval(rhs(eval(eq17, eqfcomplex)), eqp))

u(x, y, t) = -12*(-b^4*yp^2+(a*yp+xp)^2*b^2-3*a)*b^2/(b^4*yp^2+(a*yp+xp)^2*b^2+3*a)^2

Alternatively,

> 

eqfp := dchange(eqp, eqfcomplex, [xp, yp], params = [a, b], simplify); eq17p := dchange(eqp, eq17, [xp, yp], params = [a, b], simplify); eqt := simplify(eval(eq17p, eqfp))

f(xp, yp, t) = (b^4*yp^2+(a*yp+xp)^2*b^2+3*a)/b^2

u(xp, yp, t) = (6*(diff(diff(f(xp, yp, t), xp), xp))*f(xp, yp, t)-6*(diff(f(xp, yp, t), xp))^2)/f(xp, yp, t)^2

u(xp, yp, t) = -12*(-b^4*yp^2+(a*yp+xp)^2*b^2-3*a)*b^2/(b^4*yp^2+(a*yp+xp)^2*b^2+3*a)^2

Eq 17 is the parametric equation for a line.
velocities are derivatives wrt t. (after eq 19)

Intercept is at the origin (x=0,y=0 at t=0), and slope is dy/dx (Eq 19).

> 

vx, vy := diff(eval(x, eqp), t), diff(eval(y, eqp), t); dydx := simplify(vy/vx)

-5*a^2-5*b^2, 10*a

-2*a/(a^2+b^2)

> 

``

NULL

Download Y-pde-line.mw

4 5 6 7 8 9 10 Last Page 6 of 81