Mariusz Iwaniuk

1476 Reputation

14 Badges

9 years, 308 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by Mariusz Iwaniuk

@digerdiga 

I updated my answer.Its the same answer like done by  the user Carl Love.Probably a better simplification can not be done.

Ei_ver2.mw

@digerdiga 

restart;

assume(R::real, z::real); is(-R+sqrt(R^2+z^2) < 0);

#false

@digerdiga 

`assuming`([convert(signum(0, R-sqrt(R^2+z^2), 1), piecewise)], [R-sqrt(R^2+z^2) < 0]);

# -1

or:

`assuming`([signum(0, R-sqrt(R^2+z^2), 1)], [R-sqrt(R^2+z^2) < 0])

#-1

Maple have some limitation about regarding the resolution of the case and falsehood.

restart:

assume(R > 0, z > 0, R in real, z in real);

coulditbe(R-sqrt(R^2+z^2) < 0);#FAIL

is(R-sqrt(R^2+z^2) < 0);#FAIL

@one pound 

From Wikipedia, using property: Ei(1,z) = - Ei(-z) and faster converging series by Ramanujan:

For example, for z = 10,and for  40 terms  an answer correct to 16 significant figures.

series.mw

@kuwait1 

A numeric check gives a big numbers in domain x=(-5,5) and z=(-5,5).Integral is divergent.

r := 2.5;
Q := proc (x, z) options operator, arrow; evalf(Int(1/(x*r^2*cos(x-y)+z*r^2*sin(z-y))^5, y = 0 .. Pi, method = _Gquad, epsilon = 1/100000)) end proc;
seq(seq(lprint(Q(x, z)), x = -5 .. 5, .25), z = -5 .. 5, .25);

-0.3091653839e55
0.7846201518e57
0.2353821795e58
-0.7841865263e57
-0.7756899663e57
0.1123370081e58
...

@ssllys You are welcome :)

@mmcdara . Thanks for checking :)

@pr0t3ctabl3 

All  my code in answer works on Maple 2017.3 and 2018. I don't have earlier version of Maple, so I will not check it.

@pr0t3ctabl3 

I edited the answer.It should work now.

All code was exectuted in Maple 2018.

If you discovered a strange result produced by Maple,
it would be useful to exhibit the simplest version where
this abnormal result is still present.Can you add actual code so people can copy-paste easily?.

@Gabriel samaila 


 

restart

H[a] := 5; P[m] := .8; V[0] := 1; N[b] := .1; P[r] := 10; N[t] := .1; B[r] := 1; S[c] := 1

5

 

.8

 

1

 

.1

 

10

 

.1

 

1

 

1

(1)

dsys := {diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-H[a]*x[3](y)-V[0]*x[3](y)-B[r]*x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -H[a]*P[m]*x[3](y)-V[0]*P[m]*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -N[t]*x[7](y)^2-V[0]*P[r]*x[7](y)-N[b]*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -V[0]*S[c]*x[9](y)+N[t]^2*x[7](y)^2/N[b]+N[t]*P[r]*x[7](y)/N[b]+N[t]*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, x[9](0) = -.68}

{diff(x[1](y), y) = 1, diff(x[2](y), y) = x[3](y), diff(x[3](y), y) = -x[6](y)-6*x[3](y)-x[8](y), diff(x[4](y), y) = x[5](y), diff(x[5](y), y) = -4.0*x[3](y)-.8*x[5](y), diff(x[6](y), y) = x[7](y), diff(x[7](y), y) = -.1*x[7](y)^2-10*x[7](y)-.1*x[7](y)*x[9](y), diff(x[8](y), y) = x[9](y), diff(x[9](y), y) = -x[9](y)+.1000000000*x[7](y)^2+10.00000000*x[7](y)+.1*x[7](y)*x[9](y), x[1](0) = 0, x[2](0) = 0, x[3](0) = 0, x[4](0) = 0, x[5](0) = 0, x[6](0) = 0, x[7](0) = -1, x[8](0) = 1, x[9](0) = -.68}

(2)

sol := dsolve(dsys, numeric, method = classical[rk4], stepsize = .1)

proc (x_classical) 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_classical) else _xout := evalf(x_classical) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _yn, _yl, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = .1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _octl := Array(1..9, {(1) = 9, (2) = 0., (3) = 0., (4) = .1, (5) = 50000, (6) = 0, (7) = 4, (8) = 1, (9) = 0}); _n := trunc(_ctl[1]); _yini := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0., (7) = 0., (8) = -1., (9) = 1.}); _y0 := Array(0..9, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0., (7) = 0., (8) = -1., (9) = 1.}); _yl := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _yn := Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}); _fcn := proc (N, X, Y, YP) option `[Y[1] = x[1](y), Y[2] = x[2](y), Y[3] = x[3](y), Y[4] = x[4](y), Y[5] = x[5](y), Y[6] = x[6](y), Y[7] = x[7](y), Y[8] = x[8](y), Y[9] = x[9](y)]`; YP[1] := 1; YP[3] := -Y[6]-6*Y[3]-Y[8]; YP[5] := -4.0*Y[3]-.8*Y[5]; YP[7] := -.1*Y[7]^2-10*Y[7]-.1*Y[7]*Y[9]; YP[9] := -Y[9]+.1000000000*Y[7]^2+10.00000000*Y[7]+.1*Y[7]*Y[9]; YP[2] := Y[3]; YP[4] := Y[5]; YP[6] := Y[7]; YP[8] := Y[9]; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "classical" elif _xout = "numfun" then return round(_ctl[6]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[3]-_y0[0] <> 0. then _xout := _ctl[3] elif _ctl[2]-_y0[0] <> 0. then _xout := _ctl[2] else error "no information is available on last computed point" end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) 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); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n), op(_pars)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[3] = 0. then [_ctl[3], seq(_yn[_i], _i = 1 .. _n)] elif not _reinit and _xout-_ctl[2] = 0. then [_ctl[2], seq(_yl[_i], _i = 1 .. _n)] else if _ctl[2]-_y0[0] = 0. or sign(_ctl[2]-_y0[0]) <> sign(_xout-_ctl[2]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/classical`(_fcn, var(_ctl), _y0[0], var(_yl), var(_yn), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), var(Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})))) catch: userinfo(2, `dsolve/debug`, print(`Exception in classical:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) else _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end if end try else try _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..9, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0}), Array(1..3, 1..9, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0})) catch: _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end try end if; if type(_errcd, 'list') or _errcd < 0 then userinfo(2, {dsolve, `dsolve/classical`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/classical`}, ` t =`, _ctl[2]); userinfo(2, {dsolve, `dsolve/classical`}, ` y =`, _yl[1]); for _i from 2 to _n do userinfo(2, {dsolve, `dsolve/classical`}, `	 `, _yl[_i]) end do; if type(_errcd, 'list') then error op(_errcd) elif _errcd+1. = 0. then Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_ctl[3]) else Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, unknown error code returned from classical %2", evalf[8](_ctl[3]), trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_yn[_i], _i = 1 .. _n)] end if end proc, (2) = Array(0..0, {}), (3) = [y, x[1](y), x[2](y), x[3](y), x[4](y), x[5](y), x[6](y), x[7](y), x[8](y), x[9](y)], (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_classical, ["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_classical, '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_classical, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_classical, '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_classical), 'string') = rhs(x_classical); 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_classical), 'string') = rhs(x_classical)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_classical) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_classical) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

(3)

plots:-odeplot(sol, [[y, x[8](y)]], y = 0 .. 2, color = ["Red"])

 

plots:-odeplot(sol, [seq([y, x[j](y)], j = 1 .. 8)], y = 0 .. 2, view = [0 .. 2, -1 .. 1], color = ["Red", "Blue", "Black", "Green", "Yellow", "Pink", "Nautical GrayViolet", "Niagara DarkOrchid"])

 

``


 

Download testODE2.mw

 

Your_Way.mw

Thanks.@Carl Love 

First 15 16 17 18 19 20 21 Last Page 17 of 30