wg&ddlmZddlmZGddeZdZd dZee_ed k(rddlZejyy) )bisect)xrangec eZdZy) ODEMethodsN)__name__ __module__ __qualname__Y/home/mcse/projects/flask/flask-venv/lib/python3.12/site-packages/mpmath/calculus/odes.pyrrsr rc |jd| x}}t|}|g} |g} |} |} |j} | d|zz|_t|D]_}|| | }t t| Dcgc]}| ||||zz} }| |z } | j | | j | at|Dcgc]}g}}t|dzD]}dg|z}d|dzz}d}t|dzD]<}t|D]}||xx|| ||zz cc<|||z dzz| z}|dz }>|| z|j |z }t|D]$}|||z||<||j ||& | |_|j}|D]3}|ds t||j|t|dz |}5|dz}|||zfScc}wcc}w#| |_wxYw)Nrr) ldexplenprecrangerappendfaconeminnthrootabs)ctxderivsx0y0tol_precnhtoldimxsysxyorigifxydserjsbkscaleradiustss r ode_taylorr4sBiiH9%%A b'C B B A A 88D1:q AA,C(.s1v711aAh7A7 FA IIaL IIaL  !:&ar&&qs $ACAQAA1Q3Z s)AaDA1aL(D)!A#a%[qb)Q   Gcggaj(E3Z $te|!A ad# $ $WWFB b6SRV_a!@AFB aKF 6 >98 's*:G&2G8G&? G!B x_1` is also fast. **Examples of first-order ODEs** We will solve the standard test problem `y'(x) = y(x), y(0) = 1` which has explicit solution `y(x) = \exp(x)`:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> f = odefun(lambda x, y: y, 0, 1) >>> for x in [0, 1, 2.5]: ... print((f(x), exp(x))) ... (1.0, 1.0) (2.71828182845905, 2.71828182845905) (12.1824939607035, 12.1824939607035) The solution with high precision:: >>> mp.dps = 50 >>> f = odefun(lambda x, y: y, 0, 1) >>> f(1) 2.7182818284590452353602874713526624977572470937 >>> exp(1) 2.7182818284590452353602874713526624977572470937 Using the more general vectorized form, the test problem can be input as (note that *f* returns a 1-element vector):: >>> mp.dps = 15 >>> f = odefun(lambda x, y: [y[0]], 0, [1]) >>> f(1) [2.71828182845905] :func:`~mpmath.odefun` can solve nonlinear ODEs, which are generally impossible (and at best difficult) to solve analytically. As an example of a nonlinear ODE, we will solve `y'(x) = x \sin(y(x))` for `y(0) = \pi/2`. An exact solution happens to be known for this problem, and is given by `y(x) = 2 \tan^{-1}\left(\exp\left(x^2/2\right)\right)`:: >>> f = odefun(lambda x, y: x*sin(y), 0, pi/2) >>> for x in [2, 5, 10]: ... print((f(x), 2*atan(exp(mpf(x)**2/2)))) ... (2.87255666284091, 2.87255666284091) (3.14158520028345, 3.14158520028345) (3.14159265358979, 3.14159265358979) If `F` is independent of `y`, an ODE can be solved using direct integration. We can therefore obtain a reference solution with :func:`~mpmath.quad`:: >>> f = lambda x: (1+x**2)/(1+x**3) >>> g = odefun(lambda x, y: f(x), pi, 0) >>> g(2*pi) 0.72128263801696 >>> quad(f, [pi, 2*pi]) 0.72128263801696 **Examples of second-order ODEs** We will solve the harmonic oscillator equation `y''(x) + y(x) = 0`. To do this, we introduce the helper functions `y_0 = y, y_1 = y_0'` whereby the original equation can be written as `y_1' + y_0' = 0`. Put together, we get the first-order, two-dimensional vector ODE .. math :: \begin{cases} y_0' = y_1 \\ y_1' = -y_0 \end{cases} To get a well-defined IVP, we need two initial values. With `y(0) = y_0(0) = 1` and `-y'(0) = y_1(0) = 0`, the problem will of course be solved by `y(x) = y_0(x) = \cos(x)` and `-y'(x) = y_1(x) = \sin(x)`. We check this:: >>> f = odefun(lambda x, y: [-y[1], y[0]], 0, [1, 0]) >>> for x in [0, 1, 2.5, 10]: ... nprint(f(x), 15) ... nprint([cos(x), sin(x)], 15) ... print("---") ... [1.0, 0.0] [1.0, 0.0] --- [0.54030230586814, 0.841470984807897] [0.54030230586814, 0.841470984807897] --- [-0.801143615546934, 0.598472144103957] [-0.801143615546934, 0.598472144103957] --- [-0.839071529076452, -0.54402111088937] [-0.839071529076452, -0.54402111088937] --- Note that we get both the sine and the cosine solutions simultaneously. **TODO** * Better automatic choice of degree and step size * Make determination of Taylor series convergence radius more robust * Allow solution for `x < x_0` * Allow solution for complex `x` * Test for difficult (ill-conditioned) problems * Implement Runge-Kutta and other algorithms r g@(Tc||dgSNrr )r&r'F_s r zodefun..s"Q!+r Fc Z|Dcgc]}j|ddd|c}Scc}w)Nr)polyval)r,ar.rs r mpolyvalzodefun..mpolyvals*145A AddGQ'555s(c8 |krtt |}|t kr |dz S d\}}} rtd||fz |||z }|}t || \}} j | j |||f||kr dSl)Nrrz$Computing Taylor series for [%f, %f]) ValueErrorrrprintr4r)r&r r,xaxbr'Frdegreer@series_boundaries series_datarverbosers r get_serieszodefun..get_seriess r6  $a ( s$% %qs# #%b/KCR<BxGHbe$AB aQ&AGC  $ $R (   R} -Bw"2&r cj|}j} _|\}}} |||z }|_ r|Dcgc]}|c}S|dS#|_wxYwcc}wr:)convertr) r&r(r,rDrEr'ykrrKr@ return_vectorworkprecs r interpolantzodefun..interpolant sy KKNxx CH$Q-KCRad#ACH "#$BRC$ $aD5L CH$sA A+ A()intlogrdpsr TypeErrorr4)rrFrrr"rGmethodrJr,rErQr;rKr@rOrHrIrrPs``` ` ` @@@@@@@@r odefunrW3sf Q'(+88B;  .C#'' " --Fxx"}H B aR6:GCRR=/K6''$   W   &T s1 B88CC__main__)NNtaylorF) r libmp.backendrobjectrr4rWrdoctesttestmodr r r r^sJ"  *XgR  zGOOr