7 years, 6 days

## Alternative...

to @C_R

 >
 >
 >
 >
 >
 >
 > restart
 > assume(-1 <= a, a < 1)
 > solve(cos(x) > a, x) assuming -1 <= a, a < 1
 > solve(cos(x) > a, x) assuming a > -1
 > solve(cos(x) > -1/2, x)
 (1)
 > # Given a value of a, f(a) returns the domain D, not necessarily connected, such # that any x ibn D verifies cos(x) < a f := proc(a) solve(cos(x) > a, x) end proc
 (2)
 > f(-1); f(0); f(1);
 (3)
 > # Procedure g belongs basicallly how procedure f does. # The only difference is the form the domain D is returned. # # Procedure g is most aimed at ploting the domains D(a) (vertical axis) # versus a. g := proc(a)   local s := solve(cos(x) > a, x):   map(u -> [op(op(1, u)), op(op(2, u))], [s]) end proc
 (4)
 > plots:-display(   seq(seq(plot([[a, g(a)[k][1]], [a, g(a)[k][2]]], color=red, thickness=7), k=1..nops(g(a))), a in [seq](-1..1, 0.01))   , labels=["a", "Domain"]   , labeldirections=[default, "vertical"] )
 >

## Maybe you should explain in simple words...

what you want to achieve in this worksheet.
Because this is not a matter of solve or dsolve: for what I understand the last two commands of your worksheet are meaningless.

Sam_mcdara.mw

## Re-parameterize the problem...

in such a way that the interval (0, +oo) is mapped onto some finite interval, for instance (0, 1).

 > restart:
 > with(plots):
 > eq1 := diff(f(x), x, x, x)+(1/2)*cos(alpha)*x*(diff(f(x), x, x))+(1/2)*sin(alpha)*f(x)*(diff(f(x), x, x)) = 0;
 > eq2 := diff(g(x), x, x)+diff(g(x), x)+(diff(g(x), x))*(diff(h(x), x))+cos(alpha)*x*(diff(g(x), x))+sin(alpha)*f(x)*g(x) = 0;
 > eq3 := diff(g(x), x, x)+diff(h(x), x, x)+1/2*(cos(alpha)*x+sin(alpha)*f(x)) = 0;
 > ics := f(0)=0, D(f)(0)=1, (D@@2)(f)(0)=a[1], g(0)=1, D(g)(0)=a[2], h(0)=1, D(h)(0)=a[3];
 > #bcs:=f(x) , g(x), h(x) tends to 0 as x tends to infinity
 (1)

ICs ONLY

 > sys := {eq1, eq2, eq3, ics}: # dsolve(sys); # It's unlikely that a formal solution does exist
 > params := convert(indets(sys, name) minus {x}, list)
 (2)
 > sol := dsolve(sys, numeric, parameters=params);
 (3)
 > # example of use: # # Step 1: instanciate parameters data := [\$1..4]: sol(parameters=data); # Step 2 (illustration): plot some quantity display(   odeplot(sol, [x, f(x)], x=0..1, color=red, legend=typeset('f(x)'))   , odeplot(sol, [x, g(x)], x=0..1, color=green, legend=typeset('g(x)'))   , odeplot(sol, [x, h(x)], x=0..1, color=blue, legend=typeset('h(x)')) )

ICs & BCs

What ICs do you want to replace by the BCs ?

A NAIVE APPROACH WHICH DOESN'T ALWAYS WORK

 > # Method 1 # # Assuming you want to replace the ICs which contain the higher derivative order, # and assuming that L is a reasonnable approximation of +oo L := 10: `ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(L)=0, g(L)=0, h(L)=0; # REMARK: # # As soon as your differential susyemen containsBCs, the use of option "parameters" # indsolve/numerics is forbidden. # Thus you have to instanciate the parameters BEFORE using dsolve. SYS := {eq1, eq2, eq3, `ICs+BCs`}: SOL := dsolve(eval(SYS, params=~data), numeric)
 (4)
 > # Plots display(   odeplot(SOL, [x, f(x)], x=0..L, color=red, legend=typeset('f(x)'))   , odeplot(SOL, [x, g(x)], x=0..L, color=green, legend=typeset('g(x)'))   , odeplot(SOL, [x, h(x)], x=0..L, color=blue, legend=typeset('h(x)')) )
 > # REMARK: # # This naive method doesn't always work L := 100: `ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(L)=0, g(L)=0, h(L)=0; SYS := {eq1, eq2, eq3, `ICs+BCs`}: SOL := dsolve(eval(SYS, params=~data), numeric, maxmesh=10^4)

A BETTER APPROACH: RE-PARAMETERIZE THE PROBLEM

 > # Method 2 # # Re-parameterize the problem by using a 1-to-1 map from (0, +oo) into (0, 1). # # Here I use the change phi := x -> 1/(1-x) in such a way that +oo is mapped to 1/ edo := eval(eqs, {f(x)=f(phi(x)), g(x)=g(phi(x)), h(x)=h(phi(x))}): convert(select(has, indets(edo, function), D), list): Changes := % =~ subsop~(-1=z, %): eval(eval(edo, Changes), {f(phi(x))=f(z), g(phi(x))=g(z), h(phi(x))=h(z)}): EDO := (numer@lhs@simplify)~(convert(eval(eval(%, phi=(x->1/(1-x))), x=z), diff));
 (5)
 > `ICs+BCs` :=  f(0)=0, D(f)(0)=1, g(0)=1, h(0)=1, f(1)=0, g(1)=0, h(1)=0; SYS := {EDO[], `ICs+BCs`}: SOL := dsolve(eval(SYS, params=~data), numeric, maxmesh=10^4)
 (6)
 > display(   odeplot(SOL, [1/(1-z), f(z)], z=0..0.999, axis[1]=[mode=log], color=red, legend=typeset('f(x)'), labels=[x, ""])   , odeplot(SOL, [1/(1-z), g(z)], z=0..0.999, axis[1]=[mode=log], color=green, legend=typeset('g(x)'))   , odeplot(SOL, [1/(1-z), h(z)], z=0..0.999, axis[1]=[mode=log], color=blue, legend=typeset('h(x)')) )
 >

## This simple thing?...

Change the title for somthong more meaningful

 > plot(   (subs(M = 1.3015, A)-subs(M = 1.2431, A))/subs(M = 1.2431, A)*100   , x = -5 .. 0   , axes = boxed   , colour = red   , tickmarks=[default, [seq(k=sprintf("%d%%", k), k in [seq](-60..0, 10))]]   , title="(A-Aref)/Aref*100%"   , gridlines=true )
 >

## They are not ODEs...

Neither differential equation 1 nor differential equation 2 are ODEs bur simply DEs.

IMO, the correct approach to handle this situation is to proceed this way (here illustrated on case 2 for simplicity):

```ode:=diff(y(x),x)^2-(1+2*x*y(x))*diff(y(x),x)+2*x*y(x) = 0;
S := [solve(ode, diff(y(x), x))];
[1, 2 y(x) x]
map(s -> dsolve(diff(y(x), x)=s), S)
[                              / 2\]
[y(x) = x + _C1, y(x) = _C1 exp\x /]
```

By chance the solutions you got are identical to those above because the two "branches" 1 and 2*x*y(x) are easy to obtain and have simple expressions.

Let's go now to case 1. Here the "branches" have extremely complex expressions and Maple, quite logically, fails in finding the general solution.
Details are given in this attached file:

As you will see in there are four "branches" along which y(x) is governed by an ode of the form

`diff(yn(x), x) = Fn(yn(x), x), n=1..4`

You can observe that each Fn(yn(x), x) is singular at x=0 (which indicates that initial conditiopns for the form y(0)=C are inappropriate, unless you use a local series expansion to fix this problem).

From a numerical point of view dsolve/numeric seems to fail too (more precisely I wans't capable to find any "good" initial condition).

The "simpler" case I treat at the end of the file, is a case I came across during my professional work.It is not that far from your case 1 (no sqrt in my "simpler" case).

## Two ways to build such random variable...

One is to build it from a Beta random variable, the other (more complex) is to define a SubjectiveBeta distribution.

Here is an incomplete example (I didn't introduce the Skewness, Kurtosis, and other statistics into the definition of SubjectiveBeta)

 > restart:
 > with(Statistics):

WAY 1:
Define B ~ Beta(alpha, beta)
And next BSD := x__min + (x__max-x__min)*B

 > B := RandomVariable('Beta' (alpha, beta))
 (1)
 > BSD := x__min + (x__max-x__min)*B
 (2)
 > # PDF of BSD f := PDF(BSD, x)
 (3)
 > # PDF of BSD in the form given in your reference eval(f) assuming x > x__min, x__max > x__min, x < x__max: subs(x__max-x__min=h, %): map(simplify, %) assuming alpha > 0, beta > 0, h > 0: simplify(%): f__BSD := subs(h=x__max-x__min, %);
 (4)
 > # A closed form of the CDF does exist contraryto what is said in your reference F := CDF(BSD, x); print(): # After accounting for assumptions: eval(F)  assuming x > x__min, x__max > x__min, x < x__max
 (5)
 > # Direct way (I don't know how to simplify this expression) 'mean(BSD)' = simplify(simplify(Mean(BSD), GAMMA), GAMMA); # alternative way (using the linearity of the "Expectation" operator 'mean(BSD)' = x__min + (x__max-x__min)*Mean(B);
 (6)
 > # Direct way (I don't know how to simplify this expression) 'variance(BSD)' =  simplify(Variance(BSD), GAMMA): # alternative way (using the linearity of the "Expectation" operator 'variance(BSD)' = (x__max-x__min)^2*Variance(B);
 (7)
 > # Direct way (fails) 'mode(BSD)' = simplify(Mode(BSD), GAMMA); # alternative way (obvious result) 'mode(BSD)' = x__min + (x__max-x__min)*Mode(B);
 (8)

WAY 2:
Define a BSD distribution

 > B   := RandomVariable('Beta' (alpha, beta)): BSD := x__min + (x__max-x__min)*B: # f and F above are used to be copied-pasted into SubjectiveBeta f   := PDF(BSD, x); F   := CDF(BSD, x) assuming x__max > x__min
 (9)
 > SubjectiveBeta := proc(A, B, P, Q)    Distribution(      PDF=unapply(        eval(          piecewise(            (x-x__min)/(x__max-x__min) < 0, 0            , (x-x__min)/(x__max-x__min) < 1            , ((x-x__min)/(x__max-x__min))^(-1+alpha)*(1-(x-x__min)/(x__max-x__min))^(-1+beta)/Beta(alpha, beta)            , 0          )/(x__max-x__min)          , [alpha=A, beta=B, x__min=P, x__max=Q]        )        , t      ),      CDF=unapply(        eval(          piecewise(            (x-x__min)/(x__max-x__min) < 0, 0            , (x-x__min)/(x__max-x__min) < 1,                ((x-x__min)/(x__max-x__min))^alpha                *                hypergeom([alpha, 1-beta], [alpha+1], (x-x__min)/(x__max-x__min))/(Beta(alpha, beta)*alpha)            , 1          )          , [alpha=A, beta=B, x__min=P, x__max=Q]        )        , t       ),       Mean = eval(                x__min+(x__max-x__min)*alpha/(alpha+beta)                , [alpha=A, beta=B, x__min=P, x__max=Q]              ),       Variance = eval(                    (x__max-x__min)^2*alpha*beta/((alpha+beta)^2*(alpha+beta+1))                    , [alpha=A, beta=B, x__min=P, x__max=Q]                  ),       Mode = eval(                x__min+(x__max-x__min) * Mode('Beta'(A, B))                , [alpha=A, beta=B, x__min=P, x__max=Q]              ),       Conditions = [A > 0, B > 0, Q > P],       RandomSample = proc(N::nonnegint)                        P +~ (Q-P) *~ Sample('Beta'(A, B), N)                      end proc    ) end proc:
 > BSD := (alpha, beta, xmin, xmax) -> RandomVariable(SubjectiveBeta(alpha, beta, xmin, xmax)): SubjBeta := BSD(alpha, beta, x__min, x__max);
 (10)
 > PDF(SubjBeta, x); CDF(SubjBeta, x); Mean(SubjBeta); Variance(SubjBeta); Mode(SubjBeta);
 (11)
 > SubjBetaNum := BSD(3, 1, -1, 4); PDF(SubjBetaNum, x); CDF(SubjBetaNum, x); Mean(SubjBetaNum); Variance(SubjBetaNum); Mode(SubjBetaNum);
 (12)
 > Histogram( Sample(SubjBetaNum, 1000) )
 >

It would have been possible to build a module instead of a procedure to define SubjectiveBeta but I don't know what is the most robust approach (I let the specialists answer this and complete my answer).

## Several errors and ill posed problems...

Here is a corrected version, but as you provide more bcs than required it's up to you to clean it up.

 >
 > kernelopts(version);
 (1)
 >
 >
 >
 >
 (2)
 >
 (3)
 >
 (4)
 >
 (5)
 > # What you want to "dsolve" is: ToDsolve := [ subs(a1,eq1), bcs, subs(a1,eq2), bcs1 ]: print~(ToDsolve):
 (6)
 > # ToDsolve contains some formal quantities: FQ := convert(indets(ToDsolve, name) minus {f, theta, x}, list)
 (7)
 > # Give them numerical values and look to what happens dsolve(eval(ToDsolve, FQ=~1), numeric)
 > # The message means that you have 3 bcs for theta as theta verifies a 2nd order ode # It's upt to you to fix that. # I just show you that once fixed thete is no longer any error: # # Attempt 1 attempt := subsop(7=NULL, ToDsolve): print~(attempt ): dsolve(eval(attempt , FQ=~1), numeric)
 (8)
 > # Attempt 2 attempt := subsop(8=NULL, ToDsolve): print~(attempt ): dsolve(eval(attempt , FQ=~1), numeric)