API tools faq
paste
Advertisement
ArXen42

EDU simple solution

ArXen42
Sep 27th, 2016
2,817
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
MapBasic 11.10 KB | None | 0 0
raw download clone embed print report
  1. posx(t) =
  2. proc (t)
  3.     local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i;
  4.     if 1 < nargs then error "invalid input: too many arguments" end if;
  5.     _EnvDSNumericSaveDigits := Digits; Digits := 15;
  6.     if _EnvInFsolve = true then
  7.         _xout := evalf[_EnvDSNumericSaveDigits](t)
  8.     else
  9.         _xout := evalf(t)
  10.     end if;
  11.     _dat := Array(1..4, {1 = proc (_xin) local _xout, _n, _y0, _yn, _yl, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par;
  12.     _xout := _xin; _ctl := Array(1..9, [4,0.,0.,.5e-2,50000,0,0,1,7]);
  13.     _octl := Array(1..9, [4,0.,0.,.5e-2,50000,0,0,1,7]);
  14.     _n := trunc(_ctl[1]); _yini := Array(0..11, [0.,sposx,svelx,sposy,svely,undefined,undefined,undefined,undefined,undefined,undefined,undefined]);
  15.     _y0 := Array(0..11, [0.,sposx,svelx,sposy,svely,undefined,undefined,undefined,undefined,undefined,undefined,undefined]);
  16.     _yl := Array(1..11, [0,0,0,0,0,0,0,0,0,0,0]); _yn := Array(1..11, [0,0,0,0,0,0,0,0,0,0,0]);
  17.     _fcn :=
  18.     proc (N, X, Y, YP)
  19.         option `[Y[1] = posx(t), Y[2] = diff(posx(t),t), Y[3] = posy(t), Y[4] = diff(posy(t),t)]`;
  20.         if (evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[10]-Y[1])-Y[2])^2+(evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[11]-Y[3])-Y[4])^2 < 0 then
  21.             YP[1] := undefined; return 0
  22.         end if;
  23.        
  24.         if (Y[10]-Y[1])^2+(Y[11]-Y[3])^2 < 0 then
  25.             YP[1] := undefined;
  26.             return 0
  27.         end if;
  28.        
  29.         if Y[5] < 0 then YP[1] := undefined; return 0 end if;
  30.        
  31.         YP[2] := Y[5]*evalf(1/((evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[10]-Y[1])-Y[2])^2+(evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[11]-Y[3])-Y[4])^2)^(1/2))*(evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[10]-Y[1])-Y[2]);
  32.         YP[4] := Y[5]*evalf(1/((evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[10]-Y[1])-Y[2])^2+(evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[11]-Y[3])-Y[4])^2)^(1/2))*(evalf(Y[5]^(1/2))*evalf(1/((Y[10]-Y[1])^2+(Y[11]-Y[3])^2)^(1/4))*(Y[11]-Y[3])-Y[4]);
  33.         YP[1] := Y[2]; YP[3] := Y[4]; 0
  34.     end proc;
  35.     _pars := [a = a, sposx = sposx, sposy = sposy, svelx = svelx, svely = svely, tgtx = tgtx, tgty = tgty];
  36.        
  37.         if not type(_xout,'numeric) then
  38.             if member(_xout,["start", "left", "right"]) then
  39.                 return _y0[0] elif _xout = "method" then
  40.                     return "classical" elif _xout = "numfun" then
  41.                         return round(_ctl[6]) elif _xout = "initial" then
  42.                             return [seq(_yini[_i],_i = 0 .. _n)] elif _xout = "parameters" then
  43.                                 return [seq(_yini[_n+_i],_i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then
  44.                                     return [seq(_yini[_i],_i = 0 .. _n)], [seq(_yini[_n+_i],_i = 1 .. nops(_pars))] elif _xout = "last" then
  45.                                         if _ctl[3]-_y0[0] <> 0. then
  46.                                             _xout := _ctl[3] elif _ctl[2]-_y0[0] <> 0. then
  47.                                                 _xout := _ctl[2]
  48.                                             else
  49.                                                 error "no information is available on last computed point" end if
  50.                                             elif _xout = "enginedata" then
  51.                                                 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 < 7 and `dsolve/numeric/checkglobals`(7,table([(1, "value")=undefined,(2, "name")=sposx,(6, "name")=tgtx,(3, "local")=sposy,(3, "name")=sposy,(1, "local")=a,(4, "value")=undefined,(6, "local")=tgtx,(1, "name")=a,(5, "name")=svely,(4, "local")=svelx,(7, "value")=undefined,(6, "value")=undefined,(2, "value")=undefined,(5, "local")=svely,(4, "name")=svelx,(7, "local")=tgty,(3, "value")=undefined,(5, "value")=undefined,(2, "local")=sposx,(7, "name")=tgty]),_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..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),var(Array(1..3, 1..4, [[0,0,0,0],[0,0,0,0],[0,0,0,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..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..3, 1..4, [[0,0,0,0],[0,0,0,0],[0,0,0,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..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..4, [0,0,0,0]),Array(1..3, 1..4, [[0,0,0,0],[0,0,0,0],[0,0,0,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`},`\t `,_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(1..5, [18446884261893717438,18446884261893717702,18446884261893717878,18446884261893709910,18446884261893710086])), 3 = [t, posx(t), diff(posx(t),t), posy(t), diff(posy(t),t)], 4 = [a = a, sposx = sposx, sposy = sposy, svelx = svelx, svely = svely, tgtx = tgtx, tgty = tgty]}); _solnproc := _dat[1]; _pars := map(rhs,_dat[4]); if not type(_xout,'numeric') then if member(t,["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(t,'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(t,["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(t,'string')); if type(_res,'list') then return _res[2] else return NULL end if elif member(t,["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(t,'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i],_i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i],_i = 1 .. nops(_pars))] end if elif type(_xout,`=`) and member(lhs(_xout),["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(t),'string') = rhs(t); if lhs(_xout) = "initial" then if type(rhs(_xout),'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 2, rhs(_xout)]) end if elif not type(rhs(_xout),'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 2, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[2] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i],_i = 1 .. nops(_pars))] else return [_res[2], 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(t),'string') = rhs(t)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(t) else _ndsol := `tools/gensym`("posx(t)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"),_Inert_EXPSEQ(ToInert(_ndsol),_Inert_VERBATIM(pointto(_dat[2][2])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol),_Inert_EXPSEQ(ToInert(t)))) end if end if; try _res := _solnproc(_xout); _res[2] catch: error end try end proc;
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!