4 years, 72 days

## A solution among other `y'`&nbs...

A solution among other
`y'`

## The simplest way to eliminate de O(...) ...

The simplest way to eliminate de O(...) term in a Taylor expansion (for instance) is to convert the result into polynom.
For instance
T := taylor(sin(t), x):
convert(T, 'polynom'):

## Hi, I think you did something unnec...

Hi,
I think you did something unnecessary complicated.
More of this generating 10000 subscribers on a single proc can be done in 11 to 12 seconds only.

 > restart:
 > with(Statistics):
 > n_draw      := 10000: prior_rate  := Sample(Uniform(0, 1), n_draw): n_trials    := 16: subscribers := CodeTools:-Usage( Vector[row](n_draw, i -> Sample(Binomial(n_trials, prior_rate[i]), 1)) ):
 memory used=0.67GiB, alloc change=492.01MiB, cpu time=11.24s, real time=10.27s, gc time=1.61s
 > Histogram(subscribers, discrete=true, view=[-1..n_trials+1, default], thickness=10); > post_rate := map(i -> if subscribers[i]=6. then prior_rate[i] end if, [\$1..n_draw]): Histogram(post_rate, view=[0..1, default]); >

But the true thing is that one can easily determine the exact distribution of post_rate:

 > restart:
 > with(Statistics):
 > n_draw      := 10000: n_trials    := 16: # One can prove subscribers is a discrete uniform random variable over {0, ..., 16} X := RandomVariable(Binomial(16, p)); P := Probability(X=K); S := int(P, p=0..1) assuming K >= 0;  # probability that "subscribers" take the value K   (1)
 > post_rate := Probability(X=6); plot(%, p=0..1)  > # If you really want to generate a sample of subscribers subscribers := CodeTools:-Usage( LinearAlgebra:-RandomVector(n_draw, generator=0..n_trials) ); Histogram(subscribers, discrete=true, view=[-1..n_trials+1, default], thickness=10):
 memory used=81.07KiB, alloc change=0 bytes, cpu time=1000.00us, real time=0ns, gc time=0ns (2)
 >

## @hajnayeb Hi again, This reply...

@hajnayeb

Hi again,

So do not take it as a final answer to your problem but as a quick survey of what can be done.

PS: I did not use a white noise but a slightly colored one : the third argument in procedure KG (choosen to 0.1) is the correlation length of the random process [should be 0 for a true white noise but this choice will not do correct results with the attached file]).

 > restart:
 > with(LinearAlgebra): with(plots): with(Statistics):
 > # Step 1: compute the "impulse response function" (IRF) defined as the solution of # the ode with a Dirac source term at time t. # To simplify I assume here the ic for the original problem are x(0)=D(x)(0)=0. # # It can be shown this IRF is the solution of an homogeneous ode complemented with # non homogeneous ic: homogeneous_edo    := m*diff(x(t),t\$2)+c*diff(x(t),t)+k*x(t)=0; non_homogeneous_ic := x(0)=0, D(x)(0)=1/m; IRF                := unapply( rhs(dsolve([homogeneous_edo, non_homogeneous_ic],x(t))), t);   (1)
 > # Step 2: Let F(t) the random excitation and let mu__F(t) its mean and R__F(t__1, t__2) it's # autocorrelation function. # Here are some results # Response X(t) of the system under the random loading F(t) X(t) = Int(F(tau)*'IRF'(t-tau), tau=-infinity..+infinity); (2)
 > # Mean value of the response mu__X(t) = Int(mu__F(tau)*'IRF'(t-tau), tau=-infinity..+infinity); mu__X(t) = Int(mu__F(t-tau)*'IRF'(tau), tau=-infinity..+infinity);  #equivalently # if F(t) is stationnary so is X(t) mu__X = mu__F*Int('IRF'(tau), tau=-infinity..+infinity);   (3)
 > # Autocorrelation function of X(t) R__X(t__1, t__2) = Int(Int(R__X(tau__1, tau__2)*'IRF'(t__1-tau__1)*'IRF'(t__2-tau__2), tau__1=-infinity..+infinity), tau__1=-infinity..+infinity); (4)
 > # Spectral density S__X(xi) # Let H(xi) the Fourier transform of IRF(tau) and S__F(xi) the spectral density of F(t) then S__X(xi) = conjugate(H(xi))*H(xi)*S__F(xi); (5)
 > # How to draw trajectories: a notional example # F(t) is a stationnary gaussian random process (look how its point realizations F are constructed), # with gaussian correlation function (look the expression of the variable K in procedure KG) # (be free to ask me all precisions you would need). # Procedure KG generates a random realization of F(t). # For each such realization there exist a trajectory X(t), solution of the ode with random # excitation KG(...). # Nb_traj=10 such trajectories are drawn KG := proc(X, Y, psi, sigma)   local NX, DX, K, mu, k, y:   NX := numelems(X);   DX := < seq(Vector[row](NX, X), k=1..NX) >^+:   K  := (sigma, psi) -> evalf( sigma^2 *~ exp~(-((DX - DX^+) /~ psi)^~2) ):   mu := add(Y) / NX;   k  := (t, sigma, psi) -> evalf( convert(sigma^2 *~ exp~(-((t -~ X ) /~ psi)^~2), Vector[row]) ):   y  := mu + k(t, sigma, psi) . (K(sigma, psi))^(-1) . convert((Y -~ mu), Vector[column]):   return y end proc:
 > NX := 10: T  := [\$1..NX]: Nb_traj := 10: for traj from 1 to Nb_traj do   F  := convert( Statistics:-Sample(Normal(0, 1), NX), list):   X := dsolve(          {            m*diff(x(t),t\$2)+c*diff(x(t),t)+k*x(t)=KG(T, F, 0.1, 1), x(0)=1, D(x)(0)=0          },          numeric, parameters=[m, c, k]        ):   X(parameters=[1, 1, 1]):   pl||traj := odeplot(X, [t, x(t)], t=0..10, color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])) end do: display(seq(pl||traj, traj=1..Nb_traj)) >

## @Rouben Rostamian  @EarlA slight si...

@Rouben Rostamian  @Earl

A slight simplification which avoids using the parallel axis theorem
Ic := int(int(((x-x__c)^2+y^2), y=-f(x)..f(x)), x=0..1);

And a funny (maybe?) different point of view

 > restart;
 > with(Statistics):
 > f := sqrt(x^(2/3) - x^2) (1)
 > # Normalization constant C := int(x, x=0..1); (2)
 > # identify f/C to the pdf of some random variable X pdf := unapply(piecewise(x<0, 0, x < 1, f/C, 0), x); X   := Distribution(PDF = pdf);  (3)
 > # Note that the position of the center of mass doesn't change if you consider only # the part of the lamina defined by y >=0. # So, in Rouben's notations x__c := Mean(X); evalf(x__c);  (4)
 > # The variance of X is equal to the inertia of the lamina around the vertical axis # defined by x = x__c I__x__c := Variance(X); evalf(%)  (5)
 > # To obtain the inertia around the center of mass of the lamina, whose position is (x__c, 0) # you need to add a term I__y__c wich corresponds to Rouben's int(int(y^2, y=-f(x)..f(x)), x=0..1). # In the statistical framework int(y^2, y=-f(x)..f(x)) is the variance V__y(x) of a random variable # uniformly distributed with support [-f(x), f(x)]. # This variance can be interpreted as a function of the random variable whose distribution is X. # Takink its mean just gives the value of I__y__c g       := unapply(f, x); V__y    := eval(Variance(Uniform(-g(x), g(x))), x=RandomVariable(X)): I__y__c := Mean(V__y);  (6)
 > I__c := I__x__c + I__y__c; evalf(%);  (7)
 >

## Hi, Customize the code below as you...

Hi,

Customize the code below as you want

 > restart:
 > interface(rtablesize=10):
 > MyErrorPlot := proc(x, y, Ey, symb, symbsize, symbcol, errwidth, errcol)   # symb     : symbol corresponding tocouples (x[i], y[i])   # symbsize : its size   # symbcol  : its color   # errwidth : width of the horizontal bars   # errcol   : color of the error bars   plots:-display(     Statistics:-ScatterPlot(x, y, symbol=symb, symbolsize=symbsize, color=symbcol),     seq(plot([[x[n], y[n]-Ey[n]], [x[n], y[n]+Ey[n]]], color=errcol), n=1..numelems(x)),     seq(plot([[x[n]-errwidth, y[n]-Ey[n]], [x[n]+errwidth, y[n]-Ey[n]]], color=errcol), n=1..numelems(x)),     seq(plot([[x[n]-errwidth, y[n]+Ey[n]], [x[n]+errwidth, y[n]+Ey[n]]], color=errcol), n=1..numelems(x))   ) end proc:
 > A := [\$1..10]: B := [seq(sqrt(i), i = 1..10)]: C := [seq(i/10, i = 1..10)]: MyErrorPlot(A, B, C, circle, 15, red, 0.2, green) >

## Hu,  Probably lengthy, crude and n...

Hu,

Probably lengthy, crude and not top programming... but it seems to work

 > restart:
 > rewrite:=proc(expr)   local terms, REW, rewritten, eq, eqloc, eqs, deq, der, expo, der0, ind, D_deg, deg:   if type(expr, `+`) then     terms := [op(expr)]:   else     terms := expr:   end if:   REW := 0:   for eq in terms do     if type(eq, `*`) then       eqs := op(eq)     else       eqs := copy(eq)     end if:     rewritten := 1:     for eqloc in [eqs] do         deq       := convert(eqloc, 'D'):       if has(deq, 'D') then         if type(deq, `^`) then           der, expo := op(deq);         else           der, expo := copy(deq), 1         end if;         der0      := op(0, der);         ind       := op(op(der0));         D_deg     := op(op(0, der0));         if numelems([D_deg]) = 1 then deg:=1 else deg := D_deg end if;         rewritten := rewritten * (f[ind-deg](x))^expo:       else         rewritten := rewritten * deq:       end if:     end do;     REW := REW + rewritten   end do:   return REW: end proc:
 > eq := f(x)*diff(f(x),x\$2)^5 + diff(f[-2](x), x)*diff(f(x), x\$3); rewrite(eq);  (1)
 >

## @Violet Just a little error in your...

@Violet

Just a little error in your R code:

library("deSolve", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")

spirs <- function(t, x, parms) {
with(as.list(c(parms, x)), {
dS=L-(1-u)*b*S*P+th*R-nu*S
dI=(1-u)*b*S*P-(nu+sg+v+m+a-ph)*I
dP= m*I-g1*P
dR=(v+a)*I-(nu+th)*R

# the order of the unknonws to return must be the same that the one in x
list( c(dS, dP, dI, dR) )
})
}

x       <- c(S=1000,P=100,I=150,R=0)
parms <- c(L=0.05,b=0.002,ph=0.01,a=0.01,m=0.01,nu=0.001,th=0.005,u=0.1, v=0.1,sg=0.01, g1=0.03)
times <- 1:60
run_d <- ode(x, time, spirs, parms, method = rkMethod("rk45f"))

PS : I've done years ago a lot of comparisons between  R deSolve and Maple dsolve[numeric] (18 and 2015). and the results, as well as computational times were very close.
nevertheless, if you want to use the lsode-lsode suite of methods, I advice you to use R instead of Maple.

## Here is what you want. (PS: the help pa...

Here is what you want.

(PS: the help pages about events are not the clearer pages in Maple, in my own opinion)

 > restart
 > ode2 := diff(varphi(t), t, t)+omega^2*sin(varphi(t)); p0 := evalf((10*(1/180))*Pi); te := 6; event1 := [[diff(varphi(t), t), halt]]; ld := dsolve([eval(ode2, omega = 2*Pi), varphi(0) = p0, (D(varphi))(0) = 0], numeric, range = 0 .. te, events = event1);     (1)
 > # This first step is necessary. # Use any value of T larger that the supposed time that fired the event T :=10: ld(T): # now catch the time when the event fired FiredTime := op( ld( eventfired= ) ); (2)
 > plots:-odeplot(ld, [t, D(varphi)(t)], t=0..FiredTime) >

## Is it this that you want? (the example ...

Is it this that you want?

(the example comes from the help pages)

 > restart:
 > with(CurveFitting):
 > S := Spline([[0,0],[1,1],[2,4],[3,3]], v); (1)
 > dS := NULL: for n from 1 to nops(S)-1 do   if n::odd then     dS := dS, op(n, S)   else     dS := dS, diff(op(n, S), v)   end if: end do: dS := piecewise(dS, diff(op(-1, S), v)); (2)
 >

## Probably closer to what is donne in Matl...

Probably closer to what is donne in Matlab

 > restart:
 > test:= proc(x):
 > return x+1, x^2:
 > end proc:
 >
 > A, B := test(y): A, B (1)
 > # unassign A and B A:='A': B:='B': A, B; ``, B := test(y): A, B  (2)
 > # unassign A and B A:='A': B:='B': A, B; A, `` := test(y): A, B  (3)
 >

## Is this what you want?  ...

Is this what you want?
(assumming derivatives with respect to x and y do commute)

 > restart:
 > alias(V=V(x, y)): alias(alpha=alpha(x,y)): alias(beta=beta(x,y)):
 > cond_1 := diff(V, x)=V*alpha: cond_2 := diff(V, y)=V*beta:
 > eval(diff(V, x\$2), cond_1); eval(diff(V, x\$2), cond_2);  (1)
 > U := copy(V): for k from 1 to 2 do   U := eval(diff(U, y), cond_2); end do;  (2)
 > der := proc(r, s)   local fx, fy:   fx := proc(W, r)     local U, k:     U := copy(W);     for k from 1 to r do       U := eval(diff(U, y), cond_2);     end do;     return U:   end proc:   fy := proc(W, s)     local U, k:     U := copy(W);     for k from 1 to s do       U := eval(diff(U, y), cond_2);     end do;     return U:   end proc:   if r=0 then     if s=0 then       V     else       Diff(V, y\$s) = fy(V, s)     end if:   else     if s=0 then       Diff(V, y\$r) = fx(V, r)     else       Diff(Diff(V, x\$r), y\$s) = fx(fy(V, s), r); # or  fy(fx(V, r), s) assuming derivatives commute     end if:   end if: end proc (3)
 > printlevel := 2: for r from 0 to 2 do   for s from 0 to 2 do     der(r, s)   end do: end do;         (4)
 >

## Your function is extremely simple for it...

Your function is extremely simple for it can be rewritten this way

ff = C + f (xx) + ... f (xx[N])

where C is a constant and f(x) = x^2 - 10 * cos(2*Pi*x)  (A and B are two constants)
Note that the hessien of ff iis thus a diagonal matrix

Each component of the gradient of ff is defined by G := u -> 2*u + 20*Pi * sin(2*Pi*u) .
Thus each component of the gradient of ff  is null for exactly the same values of xx[n].
In the range xx[n]=-5..5 (for each n) the function u -> 2*u + 20*Pi * sin(2*Pi*u)  has 21 zeroes.
To find their values use this pice of code
Z := NULL:
z := -5:
k := 1;
while k <= 21 do
z := RootFinding:-NextZero(u->2*u+20*Pi*sin(2*Pi*u), z);
Z := Z, z:
k := k+1:
end do:
Z;

The 2nd derivative of u -> 2*u - 10 * cos(2*Pi*u) is  the function H := u -> 2 + 20*Pi * sin(2*Pi*u)
Note that the hessien of ff is a diagonal matrix whose diagonal elements are just H(xx[n])
The values of H(u) for each valus of Z are calculated by
Hval := map(u -> 2+20*Pi*sin(2*Pi*u), [Z])

Among these 21 values, 13 of them are poitive :
numelems(select[flatten](`>`, Hval, 0))

Then ff has a minimum if all the components (1000) of its hessian are negative at the points where the gradient of ff is null, and has a maximum if all the components of its hessian are positve at the points where the gradient of ff is null,
As 13 values of Hval are positive, one might expect 13^1000 minima (and 8^1000 maxima) for ff.
The number of saddle points is 21^1000 - 13^1000 - 8^1000

The largest value of ff at a point where its gradient is 0, is C + 1000 * V where V is the highest value of 2*u - 10 * cos(2*Pi*u) (u in [Z]).
This value V is about 18.94... and is obtained for z = 4.52299.. (twentieth zero).
It appears that H(4.52299..) is negative.
Thus ff reaches a maximum if xx = ... = xx = 4.52299..  and then takes the value
C + 18.94.. * 1000 (which is the highest value of ff within [-5, 5]^1000

If I'm not mistaken

## Here is an example (data XY come from th...

Here is an example (data XY come from the CurveFitting:-Spline help page)
The trap would be to search the extrema of each pice of the global function without accounting for the interval this piece is defined on. Another trap would be to search for extrema outsise the interval defined by the data.

 > restart:
 > with(CurveFitting):
 > XY := [[0,0],[1,1],[2,4],[3,3]]; S := Spline([[0,0],[1,1],[2,4],[3,3]], v);  (1)
 > # ranges bounds := min(op~(1, XY)), seq(rhs(op(k, S)), k in [seq(1..nops(S)-1, 2)]), max(op~(1, XY)); ranges := seq(bounds[k]..bounds[k+1], k=1..numelems([bounds])-1)  (2)
 > # 1st derivative of the S dS := seq(diff(op(k, S), v), k in [seq(2..nops(S), 2)]), diff(op(-1, S), v); (3)
 > # extrema sol := seq([fsolve(dS[k], v=ranges[k])], k=1..numelems([dS])) (4)
 > # reconstruction as a piecewise solution piecewise(seq(op([v >= op(1, ranges[k]) and v < op(2, ranges[k]), sol[k]]), k=1..numelems([sol]))) (5)
 >

## @JAMET  @KitonumHi JAMET, I do...

Hi JAMET,

I don't know if the solution Kitonum has provided suits you (I understood you were looking for something more formal?).
If I'm mistaken forget what comes below, otherwise this worksheet might interest you.
From a line L defined by its cartesian equation it contains the step by step procedure to:

1. construct the equation of L in the complex plane
2. construct the equation of the line orthogonal to L and which passes through (0, 0)
3. construct the expression of the intersection of these two lines

All of this is done formally, two examples follow with drawings

 > restart:
 > alias(conj = conjugate):
 > assumptions := (a::real, b::real, c::real, x::real, y::real) (1)
 > # Cartesian equation of line L L := a*x + b*y - c (2)
 > # let Z a complex number rel_1 := Z = x + I*y (3)
 > # let Z bar its conjugate rel_2 := `#mrow(mover(mo(Z),mo("̅")))` = eval(conj(Z), rel_1) assuming assumptions (4)
 > # explain x and y in terms of Z and Z bar interface(warnlevel=0): rel_3 := op( solve([rel_1, rel_2], [x, y]) ) (5)
 > # rewrite L by using the above equalities eval(L, rel_3); (6)
 > # take the numerator of this expression numer(%); (7)
 > # collect the result according to Z and Z bar `#mrow(mo("ℒ"))`  := collect(%, [Z, `#mrow(mover(mo(Z),mo("̅")))`]) (8)
 > # define the complex number v from a and b rel_4 := v = a+I*b (9)
 > # let v bar its conjugate rel_5 := `#mrow(mover(mo(v),mo("̅")))` = eval(conj(v), rel_4)  assuming assumptions; (10)
 > # reverse rel_4 and rel_5 rel_6 := (rhs=~lhs)~([rel_4, rel_5]) (11)
 > # here is the equation of L in the complex plane `#mrow(mo("ℒ"))` := eval(`#mrow(mo("ℒ"))`, rel_6) (12)
 > # find the complex equation of the line orthogonal to L and passing through the origin [0, 0] # # firstly define w as a rotation by Pi/2 of the vector whose affix is [a, b] in the complex plane rel_7 := w = v*I (13)
 > # let w bar its conjugate rel_8 := `#mrow(mover(mo(w),mo("̅")))` = eval(conj(w), rel_7); rel_8 := eval(rel_8, conj(v) = `#mrow(mover(mo(v),mo("̅")))`);  (14)
 > # here is the equation of the line orthogonal to L, and which passes through the origin `#mrow(mo("ℋ"))` := eval(                                 `#mrow(mo("ℒ"))`,                                 [                                   v=w,                                   `#mrow(mover(mo(v),mo("̅")))` = `#mrow(mover(mo(w),mo("̅")))`,                                   c=0                                 ]                               ) (15)
 > # where do lines L and H intersect ? # # let's form this system in Z and Z bar (brute force) sys := [ `#mrow(mo("ℒ"))` = 0, eval(`#mrow(mo("ℋ"))`, [rel_7, rel_8]) = 0 ]: print~(sys):  (16)
 > # where do lines L and H intersect ? sol := op( solve( sys, [ Z, `#mrow(mover(mo(Z),mo("̅")))` ] ) ) (17)
 > # use the expression of v to explicit the solution (obviously Z and Z bar are conjugate) sol := eval(sol, [rel_4, rel_5]); (18)
 > # transform the solution into a prettier form sol := lhs~(sol) =~ numer~(rhs~(sol)) *~ conj~(denom~(rhs~(sol)))                     /~                     expand~(denom~(rhs~(sol)) *~ conj~(denom~(rhs~(sol))) )  assuming assumptions (19)
 > # this complex number -represents the affix of the intersection point of the two lines `#mrow(mo("ℤ"))` := simplify(eval(Z, sol));  # asuming a^2+b^2 <> 0 (20)
 > # example data := [a=3, b=-2, c=1, x=X, y=Y]: with(plots): f := eval(a*X + b*Y - c, data); `#mrow(msub(mo(f),mo("⊥")))` := simplify(                                         eval(                                           eval(                                             eval(`#mrow(mo("ℋ"))`, [rel_1, rel_2, rel_7, rel_8]),                                             [rel_4, rel_5]                                           ),                                           data                                         )                                       ); __Z := eval(`#mrow(mo("ℤ"))`, data); display(   implicitplot(f, X=-1..1, Y=-1..1, color=blue),   implicitplot(`#mrow(msub(mo(f),mo("⊥")))`, X=-1..1, Y=-1..1, color=red),   pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),   gridlines=true,   scaling=constrained );    > # just another example data := [a=0, b=-3, c=-2, x=X, y=Y]: with(plots): f := eval(a*X + b*Y - c, data); `#mrow(msub(mo(f),mo("⊥")))` := simplify(                                         eval(                                           eval(                                             eval(`#mrow(mo("ℋ"))`, [rel_1, rel_2, rel_7, rel_8]),                                             [rel_4, rel_5]                                           ),                                           data                                         )                                       ); __Z := eval(`#mrow(mo("ℤ"))`, data); display(   implicitplot(f, X=-1..1, Y=-1..1, color=blue),   implicitplot(`#mrow(msub(mo(f),mo("⊥")))`, X=-1..1, Y=-1..1, color=red),   pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),   gridlines=true,   scaling=constrained );    >

 1 2 3 4 5 6 7 Last Page 1 of 13