restart

alias(theta=theta(x))

# Equation b5 with f set to 1

b5 := {1.116150082*(diff(theta, x, x))+3.100000000*(diff(theta, x)) = 0, theta(0) = 1, theta(10) = 0};
dsolve(%, numeric)

Error, (in dsolve/numeric/type_check) insufficient initial/boundary value information for procedure defined problem

alias(theta=theta)

b5 := {1.116150082*(diff(theta(x), x, x))+3.100000000*(diff(theta(x), x)) = 0, theta(0) = 1, theta(10) = 0};
dsolve(%, numeric)

alias(f=f(x), theta=theta(x));

alias(theta=theta, f=f);
b5;
bcs1; # wrong expression

# Redefine bcs1
bcs1 := theta(0) = 1, theta(10) = 0;


# 1st try (fails because f(x) is undefined)
dsolve({b5, bcs1}, numeric);

# 2nd try: give f(x) an explicit expression, for instance f(x)=x

{eval(b5, f=(x-> x)), bcs1};
dsolve(%, numeric)

Error, (in dsolve/numeric/bvp/convertsys) the ODE system does not contain derivatives of the unknown function f

# If f(x) is the solution from another ode, for instance
#
# Method 1


bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(10) = 1;
c1  := dsolve({b1, bcs}, numeric);

__f := proc(__x);
   if not type(evalf(__x),'numeric') then
      'procname'(__x);
   else
      eval(f(x), c1(__x));
   end if;
end proc:

# usage

c1(2);
'f(2)' = __f(2);
plot(__f(x), x=0..10, title="f as the solution of b1");

# Represent __f(x) by a spline function (the higher the number of discretizationpoints
# the higher the cvalue of maxmesh).

F := unapply( CurveFitting:-Spline([seq([__x, __f(__x)], __x in [seq](0..10, 0.1))], x, degree=3), x):
{eval(b5, f=(x-> F(x))), bcs1}:
c5 := dsolve(%, numeric, maxmesh=1024, abserr=1e-4):

odeplot(c5, [x, theta(x)], 0 .. 7, colour = blue, title="solution of b5 with f solution of b1")

# If f(x) is the solution from another ode, for instance
#
# Method 2: use optio known 

c5 := dsolve({eval(b5, f=__f), bcs1}, numeric, known=__f);
odeplot(c5, [x, theta(x)], 0 .. 7, colour = blue, title="solution of b5 with f solution of b1")

restart;

kernelopts(version):

de := diff(u(y), y, y) + 4*beta*diff(u(y), y)^(4 - 1)*diff(u(y), y, y)+3=0;
de := isolate(de, diff(u(y), y, y)):

de := unapply(de, beta);

desol := proc(x, k, beta, N, e)
  local sigma := 1+x^2/2:
  dsolve({de(beta), u(-sigma)=1, u(sigma)=k}, numeric, output=listprocedure, maxmesh=N, abserr=e):
end proc:

# Data for which the Newton converges

X   := 0:
r   := -(1+X^2/2)..(1+X^2/2);
sol := desol(X, 3, 3/2, 10^4, 1e-6);
f1  := eval(u(y), sol):
f2  := eval(diff(u(y), y), sol):
plots:-dualaxisplot(plot(f1(y), y=r, color=blue), plot(f2(y), y=r, color=red))

# Another data for which the Newton still converges


X   := 0:
r   := -(1+X^2/2)..(1+X^2/2);
sol := desol(X, -5, 1/2, 10^4, 1e-6);
f1  := eval(u(y), sol):
f2  := eval(diff(u(y), y), sol):
plots:-dualaxisplot(plot(f1(y), y=r, color=blue), plot(f2(y), y=r, color=red))

# And a non convergent case

X   := 0:
r   := -(1+X^2/2)..(1+X^2/2);
sol := desol(X, 1, 1/2, 10^4, 1e-6);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

DE := Diff(u(y), y, y) + 4*beta*Diff(u(y), y)^(4 - 1)*Diff(u(y), y, y)+3=0;

DE := Diff(u(y), y, y) + beta*Diff(Diff(u(y), y)^4, y)+3=0;

DE := Diff(u(y), y) + beta*(Diff(u(y), y)^4)+3*y=C;

#------------------------------------------------------------
#
# Procedures
Xsol := proc(__beta, __k, X, i)
  local data, sols, usol, csol, T, Left_bc, Right_bc, check_de:
  data := [beta=__beta, k=__k, L=1+X^2/2];
  eval({diff(u(y), y)=ade[i], u(-L)=1}, data), data
end proc:


Cstart := proc()
  local c := 0;
  local StartNotFound := true:
  while StartNotFound do
    try
      sol(parameters=[c]);
      sol(eval(L, data));
      StartNotFound := false:
    catch:
      c := c+0.1
    end try:
  end do:
  c
end proc:


FindC := proc(c)
  sol(parameters=[c]);
  eval(u(y), sol(eval(L, data)))
end proc:


f := proc(__y)
  eval(u(y), sol(__y))
end proc:

#------------------------------------------------------------

solve(DE, Diff(u(y), y));
ade := [allvalues(%)]:

# Select the 3rd solution delivered by allvalues(...)

sys, data := Xsol(1/2, 4, 0, 3):
sol := dsolve(sys, numeric, parameters=[C]);

# Find the smaller value oc C such that the solution might be
# constructed up to L.

c0 := Cstart();

# Check if c0 is admissible.

sol(parameters=[c0]):
eval(u(y), sol(eval(L, data)));
if % < eval(k, data) then
  admissible := true
else
  admissible := false
end if;


if admissible then
  # From C=c0 find the value ck of C such that y(L)=k.

  ck := fsolve('FindC(c)'=eval(k, data), c=c0..10);

  # Instanciate the solution with C=ck

  sol(parameters=[ck]):

  # Plot it.

  plot('f(__y)', __y=eval(-L..L, data))
end if;

# For k=2, and X=0 (==> 1+X^2/2 = 1) thesmallest value c9 of C
# such that the solution can be computed up to y=1 is c0=2.4.
# But y(1) = 3.12 being already larger then 2, there is no value
# of C such that y(1) = 2.

sys, data := Xsol(1/2, 2, 0, 3):
sol := dsolve(sys, numeric, parameters=[C]);

# Find the smaller value oc C such that the solution might be
# constructed up to L.

c0 := Cstart();

# Check if c0 is admissible.

sol(parameters=[c0]):
eval(u(y), sol(eval(L, data)));
if % < eval(k, data) then
  admissible := true
else
  admissible := false
end if;


if admissible then
  # From C=c0 find the value ck of C such that y(L)=k.

  ck := fsolve('FindC(c)'=eval(k, data), c=c0..10);

  # Instanciate the solution with C=ck

  sol(parameters=[ck]):

  # Plot it.

  plot('f(__y)', __y=eval(-L..L, data))
end if;

restart:

with(Statistics):

A__1 := RandomVariable(Normal(0, sigma__1));
A__2 := RandomVariable(Normal(0, sigma__2));

RV := sigma__1 * rho__12 * A__1 + sigma__2 * sqrt(1-rho__12^2) * A__2 + nu__0

# Observe that RV doesn't inherit of the type "RandomVariable"

type~([A__1, A__2, RV], RandomVariable);

# You don't whant this

expand(RV*(sigma__2 * rho__12 + nu__0));

# but this

Z := 'RV'*(sigma__2 * rho__12 + nu__0);
Z;

Z+1;

f(Z);

restart

de := diff(u__0(y), y, y) + beta*diff(u__1(y), y, y) + beta^2*diff(u__2(y), y, y) + n*beta*(diff(u__0(y), y) + beta*diff(u__1(y), y) + beta^2*diff(u__2(y), y))^(n - 1)*(diff(u__0(y), y, y) + beta*diff(u__1(y), y, y) + beta^2*diff(u__2(y), y, y)) - diff(p__0(x), x) - beta*diff(p__1(x), x) - beta^2*diff(p__2(x), x)

# This procedure returns a table whose indices are the powers of beta
# and the entries their coefficients when n is set to the value defined
# by the parameter of ExtractBeta.

ExtractBeta := proc(N)
  local c := coeffs(expand(eval(de, n=N)), beta, 'b');
  table([seq(b[i]=c[i], i=1..numelems([b]))])
end proc:

T := ExtractBeta(2)

restart:

# Let A=beta and B=beta^2, and define f this way

f := (p+q*A+r*B)^n

# The multinomial coefficient of A^s*B^t (=beta^(s+2*t)) in the expansion of f
# has expression:


MultinomialCoeff := (n, s, t) -> n!/((n-s-t)!*s!*t!) * p^(n-s-t) * q^(s) * r^(t);

# The set of all couples couples (s, t) such that s+2*t has a given value u
# can be obtained using combinat:-composition.

st_set := proc(u)
  if   u >= 3 then map(x -> [x[1], x[2]/2], select((x -> x[2]::even), combinat:-composition(u, 2))):
  elif u =  2 then {[2, 0], [0, 1]}:
  elif u =  1 then {[1, 0]}:
  else      {[0, 0]}
  end if:
end proc:

# Examples

st_set(1);
st_set(2);
st_set(6);

# Let u some positive integer.
# Compute the multinomial coefficient of beta^u in the expansion of f.

MultinomialCoeffBeta := proc(n, u)
  local MC, st_set:
  MC     := (n, s, t) -> n!/((n-s-t)!*s!*t!) * p^(n-s-t) * q^(s) * r^(t);
  st_set := proc(u)
    if   u >= 3 then map(x -> [x[1], x[2]/2], select((x -> x[2]::even), combinat:-composition(u, 2))):
    elif u =  2 then {[2, 0], [0, 1]}:
    elif u =  1 then {[1, 0]}:
    else      {[0, 0]}
    end if:
  end proc:
  add(MC(n, op(x)), x in st_set(u))
end proc:

# Examples

MultinomialCoeffBeta(n, 0);
MultinomialCoeffBeta(n, 1);
MultinomialCoeffBeta(n, 2);
MultinomialCoeffBeta(n, 6);

# Consider now the following expression:

g := (p+q*A+r*B)^n * (P+Q*A+R*B);

# How to find the expression of beta^u in the expansion of g?
# Expanding g gives:

g := expand(g);

# Then the coefficient of beta^u comes from three sources:
# source 1: this is MultinomialCoeffBeta(n, u) * P
# source 2: this is MultinomialCoeffBeta(n, u-1) * Q
# source 3: this is MultinomialCoeffBeta(n, u-2) * R

# Example

P*MultinomialCoeffBeta(n, 6) + Q*MultinomialCoeffBeta(n, 5) + R*MultinomialCoeffBeta(n, 4)

simplify(%, size);

# Your expression de is made of 3 terms:

Term_1 := diff(u__0(y), y, y) + beta*diff(u__1(y), y, y) + beta^2*diff(u__2(y), y, y);
Term_2 := n*beta*(diff(u__0(y), y) + beta*diff(u__1(y), y) + beta^2*diff(u__2(y), y))^(n - 1)*(diff(u__0(y), y, y) + beta*diff(u__1(y), y, y) + beta^2*diff(u__2(y), y, y));
Term_3 :=  - diff(p__0(x), x) - beta*diff(p__1(x), x) - beta^2*diff(p__2(x), x)

# Applying the previous rules enables computing the coefficient of beta^u in Term_2:


Coeff_2 := proc(u)
  if u = 0 then
    0:
  elif u = 1 then
   n*P*MultinomialCoeffBeta(n-1, u-1):
  elif u = 2 then
   n*(P*MultinomialCoeffBeta(n-1, u-1) + Q*MultinomialCoeffBeta(n-1, u-2)):
  else
   n*(P*MultinomialCoeffBeta(n-1, u-1) + Q*MultinomialCoeffBeta(n-1, u-2) + R*MultinomialCoeffBeta(n-1, u-3)):
  end if:
end proc:

Coeff_2(0);
Coeff_2(1);
Coeff_2(2);
Coeff_2(6);

# Define the substitution rules

Equivalences := [p, q, r, P, Q, R] =~ [coeffs(op(1, op(-2, Term_2)), beta), coeffs(op(-1, Term_2), beta)]

# And apply them to compute the "true" coefficient of betaû in Term_2
# Example 1

eval(Coeff_2(0), Equivalences);
print();
eval(Coeff_2(1), Equivalences);
print();
eval(Coeff_2(2), Equivalences);

# Example 2

eval(Coeff_2(6), Equivalences)

# For term_1, obviously:

Coeff_1 := proc(u)
  if u < 3 then coeffs(Term_1, beta)[u+1]
  else     0
  end if:
end proc:

# Examples

Coeff_1(0);
Coeff_1(1);
Coeff_1(2);
Coeff_1(3);

# At last, for Term_3

Coeff_3 := proc(u)
  if u < 3 then coeffs(Term_3, beta)[u+1]
  else     0
  end if:
end proc:

# Examples

Coeff_3(0);
Coeff_3(1);
Coeff_3(2);
Coeff_3(3);

# It remains to aggregate these three partial results

Coeff_Beta := proc(u)
  Coeff_1(u) + eval(Coeff_2(u), Equivalences) + Coeff_3(u)
end proc:

# Examples

Coeff_Beta(0);
print():
Coeff_Beta(1);
print():
Coeff_Beta(2);
print():
Coeff_Beta(3);

restart

f := sqrt(a*(a+b+c))/(a*b+a*c+b*c)

g := (a*b+a*c+b*c)*((1/3)*a+(1/3)*b+(1/3)*c)^(1/3)

map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, [op(f)])

map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, [op(g)])

map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, [op(f+g)])

map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, indets(f));
map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, indets(g));
map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, indets(f+g));

h := f^(1/2);
map(t -> if patmatch(t, a::anything^n::fraction,'p') then eval(1/n, p) end if, indets(h));

EL     := [convert~(E, list)[]]:
__from := map2(op, 1, EL):
__to   := map2(op, 2, EL):

sprintf(cat(seq("""%d"",", k=1..numelems(__from))), op(__from)):
R_from := cat("From <- c(", %[1..-2], ")");
sprintf(cat(seq("""%d"",", k=1..numelems(__to))), op(__to)):
R_to := cat("To <- c(", %[1..-2], ")");

R_ft := "ft <- cbind(From, To)";
R_gr := "graph1 <- ftM2graphNEL(ft)";

MyRsource := cat(currentdir(), "/graph.R"):

writeline(MyRsource, R_from):
writeline(MyRsource, R_to):
writeline(MyRsource, R_ft):
writeline(MyRsource, R_gr):
fclose(MyRsource)
 

restart:

f := x -> (x^2+1)/(x^2-1)

# As f is summetric the centre (Omega) of the circle C is on the vertical axis.
# Let (0, h) the coordinates of Omega.
# One has f(0)=-1.
# Let R the radius of C, then R = h+1.
#
# Let's take the branch of f defined by x > 1and search the point P of
# coordinates (X, Y=f(X)) such that (X-0)^2+(Y-h)^2=(h+1)^2:


CircleEquation = X^2+(f(X)-(R-1))^2 - R^2;

# As X > 1 on gets this relation


rel := numer(expand(rhs(%)))

# Clearly only the  term in parenthesis matters.
# Set X^2 = Z

rel1 := collect(algsubs(X^2=Z, op(-1, rel)), Z)

# As we are looking for a tangence point, the discriminant of this
# equation must be equal to 0.
discrim(rel1, Z);
solve({%=0, R > 0}, R)

# Thus R = 2 and the circle C has center (0, R-1=1)

R_circle := 2:
 

# Contact point for X > 1.
# As the squareroot ithin sol is equal to 0, one has

Z_contact := 2*R_circle - 1;

# And so:

X_contact := surd(Z_contact, 2);
Y_contact := (Z_contact+1) / (Z_contact-1);

# The contact point for X < 1 has coordinates

[-X_contact, Y_contact]