mmcdara

This answers your question and not anoth...

Answered: mmcdara 7891

March 26 2023

1 2

This answers your question and not another.

MTM:-sort([x, x/3])
                            [1     ]
                            [- x, x]
                            [3     ]

Using sort and options `>`or `<` (default) is useless because we don't know ix/3 is larger or lower than x

sort([1, 1/3]);          # 1/3, 1
sort([1, 1/3], `>`);     # 1, 1/3
sort([-1, -1/3]);        # -1, -1/3
sort([-1, -1/3], `>`);   # -1/3, -1

Hi, I guess it's more a problem...

Answered: mmcdara 7891

March 26 2023

1 8

Hi,

I guess it's more a problem of visualization than a problem of storage:

V := Vector[row](6, {1 = T[1], 2 = T[2]}, storage = sparse, fill = {}); 
op(%)
            Vector[row](%id = 18446744078202550262)
6, {1 = T[1], 2 = T[2]}, datatype = anything, storage = sparse, 

  order = Fortran_order, shape = []
addressof(V[1]);
addressof(V[2]);
addressof(V[3]);
addressof(V[4]);
addressof(V[5]);
addressof(V[6]);
                      18446744078202681118
                      18446744078202681166
                      18446744073709551615
                      18446744073709551615
                      18446744073709551615
                      18446744073709551615

The adresses of all the unassigned elements of V being the same I understand (but maybe I'm wrong) that no memory is allocated to them.
By comparison

V := Vector[row](6, symbol=s); 
addressof(V[1]);
addressof(V[2]);
addressof(V[3]);
addressof(V[4]);
addressof(V[5]);
addressof(V[6]);
                      18446744078205130166
                      18446744078205130190
                      18446744078205130214
                      18446744078205130238
                      18446744078205130262
                      18446744078205130286

Two ways depending on what you want...

Answered: mmcdara 7891

March 26 2023

2 0

UPDATED

BEFORE ALL: LogNormal distribution already exists in MAPLE.

Let us suppose that LogNormal is not a member of the ListKnownDistributions list.
There are two ways to work with a LogNormal (LN) distribution:

First way is to create a LN variable Y from a Normal one (the underlying normal rv Carl refered to).
1. define X as a Normal rv
2. set Y = exp(X)
A more complex way is to define a Distribuion of LN type (see the attached file).
Note this is purely academic for this distribution already exists in Maple: but this may help you to define new distributions (a LN distribution parameterized diffferentlyf, a normal distribution parameterized by mean and variance, a Dirichlet distribution , ...), the process is exactly the same as the one described here.

At the end of the file you will see a list of all the statistics you can define for a continuous random variable. You can add them in the definition of LogNormal to build a more complete distribution.

LogNormal.mw

I don't about such build-in function...

Answered: mmcdara 7891

March 24 2023

1 0

I don't about such build-in function but you can proceed this way

>

restart:

>	NR := 10: NC := 8: M := LinearAlgebra:-RandomMatrix(NR, NC, generator=0..1)

(1)

>	ZerosInRow := i -> NC-numelems(op(M[i])[2])

(2)

>	ZerosInRow(1); ZerosInRow(2);

(3)

>

Download ZerosInRow.mw

Answer to your 230323_Problem_Idea.mw....

Answered: mmcdara 7891

March 23 2023

2 24

@MaPal93

It's 11 pm here, so I'm going to quit Mapleprimes in a few minutes.
Just the time to provide you a worksheet which implements your ides, but do not validate your last result.
Please check what I did (all the intermediate results about the expectations are correct).

There is just an ambiguous point about the last operation you mention.

>

restart:

PART 1

>	with(Statistics):

>	X := RandomVariable(Normal(mu__X, sigma__X)): Y := RandomVariable(Normal(mu__Y, sigma__Y)):

>	Correlation(X, Y) assuming sigma__X > 0, sigma__Y > 0

(1)

>	# To correlate X and Y use a third variable Z defined this way: Z := aX + bY

(2)

>	# We want these two relations to be verified: rel_1 := Correlation(X, Z) = rho assuming sigma__x::positive, sigma__Y::positive; rel_2 := Variance(Z) = sigma__Y^2 assuming sigma__X::positive, sigma__Y::positive;

(a*mu__X^2+a*sigma__X^2+b*mu__X*mu__Y-mu__X*(a*mu__X+b*mu__Y))/(sigma__X*(a^2*sigma__X^2+b^2*sigma__Y^2)^(1/2)) = rho

(3)

>	# Solve these two relations with respect to a and b. solve([rel_1, rel_2], [a, b]); av := allvalues(%);

[[a = -rho*sigma__Y/sigma__X, b = RootOf(_Z^2+rho^2-1)], [a = rho*sigma__Y/sigma__X, b = RootOf(_Z^2+rho^2-1)]]

[[a = -rho*sigma__Y/sigma__X, b = (-rho^2+1)^(1/2)], [a = rho*sigma__Y/sigma__X, b = (-rho^2+1)^(1/2)]], [[a = -rho*sigma__Y/sigma__X, b = -(-rho^2+1)^(1/2)], [a = rho*sigma__Y/sigma__X, b = -(-rho^2+1)^(1/2)]]

(4)

>	# Let's take the first solution (this is an arbitrary decision). ab := av[1][1]

[a = -rho*sigma__Y/sigma__X, b = (-rho^2+1)^(1/2)]

(5)

>	# Finally copy Z into Y if you want to keep working with X and Y. # Don't forget to check the result. Y := copy(Z): eval([Correlation(X, Y), Variance(Y)], ab) assuming sigma__X::positive, sigma__Y::positive: simplify(%) assuming sigma__Y::positive;

(6)

APPLICATION TO YOUR PROBLEM

>	nu__jk := RandomVariable(Normal(q__0jk, sqrt(Sigma__0jk))): nu__ki := RandomVariable(Normal(q__0ki, sqrt(Sigma__0ki))):

>	# Correlate nu_jk and nu_ki (correlation coeff = rho__XX__23) XX := eval(ab, [sigma__X=sqrt(Sigma__0jk), sigma__Y=sqrt(Sigma__0ki), rho=rho__XX__23]);

[a = -rho__XX__23*Sigma__0ki^(1/2)/Sigma__0jk^(1/2), b = (-rho__XX__23^2+1)^(1/2)]

(7)

>	nu__ki := eval(anu__jk + bnu__ki, XX)

-_R7*rho__XX__23*Sigma__0ki^(1/2)/Sigma__0jk^(1/2)+_R8*(-rho__XX__23^2+1)^(1/2)

(8)

>	# WATCH OUT!!! # The expectation of nu_ki is no longer q__0ki: expand(Mean(nu__ki)); # thus we must shift nu__ki nu__ki := nu__ki - expand(Mean(nu__ki)) + q__0ki:

q__0ki*(-rho__XX__23^2+1)^(1/2)-q__0jk*rho__XX__23*Sigma__0ki^(1/2)/Sigma__0jk^(1/2)

(9)

>	# Check eval([Correlation(nu__jk, nu__ki), Variance(nu__ki), Mean(nu__ki)]) assuming Sigma__0jk::positive, Sigma__0ki::positive: simplify(%) assuming sigma__Y::positive;

(10)

>	# It's probably interesting to define some assumptions to use hereafter assumptions := Sigma__0jk::positive, Sigma__0ki::positive:

>	u__jk := RandomVariable(Normal(0, sigma__ujk)): u__ki := RandomVariable(Normal(0, sigma__uki)):

>	Variance(nu__jk + nu__ki); Mean(nu__jk^2); Mean(nu__ki^2); Mean(u__jk^2); Mean(u__ki^2);

(11)

>	allrv := ["nu__jk", "nu__ki", "u__jk", "u__ki"]: for i from 1 to 2 do rvi := parse(allrv[i]): for j from 3 to 4 do rvj := parse(allrv[j]): printf("E(%s, %s) = %a\n", allrv[i], allrv[j], Mean(rvi*rvj)): end do: end do:

E(nu__jk, u__jk) = 0
E(nu__jk, u__ki) = 0
E(nu__ki, u__jk) = 0

E(nu__ki, u__ki) = 0

>	Mean(nu__jk*nu__ki);

PART 2 (minimization problem)

>	argmin := Mean( (nu__jk(-beta__jk2lambda__jk+1)-lambda__jk(nu__kibeta__jk3+u__jk)-alpha__jk*lambda__jk-mu__jk)^2 )

(12)

>	FOC_1A := diff(argmin, mu__jk) = 0;

(13)

>	FOC_1B := diff(argmin, lambda__jk) = 0;

(14)

>	# I don't really understand your phrase # Plug FOC_1A into FOC_1B and FOC_1B can now be solved for . # I guess you want to do this (?) solve([FOC_1A, FOC_1B], [mu__jk, lambda__jk]);

(15)

>	eval(lambda__jk, %[]);

Sigma__0jk^(1/2)*(-Sigma__0ki^(1/2)*beta__jk3*rho__XX__23+Sigma__0jk^(1/2)*beta__jk2)/(-2*Sigma__0jk^(1/2)*Sigma__0ki^(1/2)*beta__jk2*beta__jk3*rho__XX__23+Sigma__0jk*beta__jk2^2+Sigma__0ki*beta__jk3^2+sigma__ujk^2)

(16)

>

# In fact "Plug FOC_1A into FOC_1B" means nothing in Maple.
# The closest way of this is eval(FOC_1B, something=something else)
# where:
# something is expression that FOC_1B contains
# something else is an expression uses as a replacementt of something.
#
# For instance

replacement := isolate(FOC_1A, mu__jk);

plug_into_FOC_1B := eval(FOC_1B, replacement);

(17)

>	solve(plug_into_FOC_1B, lambda__jk);

(Sigma__0ki^(1/2)*Sigma__0jk^(1/2)*beta__jk3*rho__XX__23-Sigma__0jk*beta__jk2)/(2*Sigma__0jk^(1/2)*Sigma__0ki^(1/2)*beta__jk2*beta__jk3*rho__XX__23-Sigma__0jk*beta__jk2^2-Sigma__0ki*beta__jk3^2-sigma__ujk^2)

(18)

>	# Which is obviously the same thing that solving a 2 equation system in mu__jk and lambda__jk

Download 230323_Problem_mmcdara.mw

combinat:-numbcomb(6, 3); ...

Answered: mmcdara 7891

March 23 2023

1 0

combinat:-numbcomb(6, 3);
                           20
combinat:-choose(6, 3);  # still tractable
   [[1, 2, 3], [1, 2, 4], [1, 2, 5], [1, 2, 6], [1, 3, 4], 

     [1, 3, 5], [1, 3, 6], [1, 4, 5], [1, 4, 6], [1, 5, 6], 

     [2, 3, 4], [2, 3, 5], [2, 3, 6], [2, 4, 5], [2, 4, 6], 

     [2, 5, 6], [3, 4, 5], [3, 4, 6], [3, 5, 6], [4, 5, 6]]
combinat:-numbcomb(100, 3);
                             161700 

CodeTools:-Usage(combinat:-choose(100, 3)):
memory used=19.90MiB, alloc change=77.36MiB, cpu time=220.00ms, real time=221.00ms, gc time=40.07ms

combinat:-choose(100, 3)[2023]
                          [1, 25, 47]

There are different ways to do this...

Answered: mmcdara 7891

March 23 2023

2 0

A few are given in the attached file

>

(1)

>

R := int(Student:-VectorCalculus:-`*`(1/e^Student:-VectorCalculus:-`*`(theta, t), int(Student:-VectorCalculus:-`*`(f(u), e^Student:-VectorCalculus:-`*`(theta, u)), u = t .. t1)), t = 0 .. t1);

(theta*ln(e)*t1-1+e^(-theta*t1))/(theta^2*ln(e)^2)

(2)

>	with(IntegrationTools):

>	# Writting < Int > instead of < int > enables using IntegrationTools, # and thus extract the integrand < c > in a simple way. # As a drawvack you must use value to evaluate the integral r := Int(1/e^(thetat)Int(f(u)e^(thetau), u = t..t1), t=0..t1); c := IntegrationTools:-GetIntegrand(r);

Int((Int(e^(-theta*u)*e^(theta*u), u = t .. t1))/e^(theta*t), t = 0 .. t1)

(Int(e^(-theta*u)*e^(theta*u), u = t .. t1))/e^(theta*t)

(3)

>

value(r);

(theta*ln(e)*t1-1+e^(-theta*t1))/(theta^2*ln(e)^2)

(4)

>	ci := IntegrationTools:-GetIntegrand(c);

(5)

>	G := unapply( Int(ci, u=0..t1) / denom(c), t)

proc (t) options operator, arrow; (Int(e^(-theta*u)*e^(theta*u), u = 0 .. t1))/e^(theta*t) end proc

(6)

>	value(G(t))

(7)

>	value(G(0))

(8)

BUT...

>	G := t -> Int(ci, u=0..t1) / denom(c); value(G(t)); value(G(0)); # You can "fix" this that way eval(value(G(t)), t=0)

proc (t) options operator, arrow; Student:-VectorCalculus:-`*`(Int(ci, u = 0 .. t1), 1/denom(c)) end proc

(9)

FINALLY you can simplify the "unapply" defiinition of G this way

>	G := unapply( value(Int(ci, u=0..t1)) / denom(c), t)

proc (t) options operator, arrow; t1/e^(theta*t) end proc

(10)

>

G(0)

(11)

>

Download M2_mmcdara.mw

A workaround...

Answered: mmcdara 7891

March 21 2023

0 4

Let's note that

u := (1 + 20230321)*x*y - (x^2 + y^2)/2

is an homogenous polynomial wrt x and y.
One can rewrite u by setting [x=lambda*cos(t), y=lambda*sin(t)] which gives

v := collect(eval(u, [x=lambda*cos(t), y=lambda*sin(t)]), lambda);
    /                         1       2   1       2\       2
    |20230322 cos(t) sin(t) - - cos(t)  - - sin(t) | lambda 
    \                         2           2        /

The maximum (maxima in reality) of v are obtained when lambda is equal to +/- infinity and when the t unction (op(1, v)) reaches its maxima:

v       := collect(eval(u, [x=lambda*cos(t), y=lambda*sin(t)]), lambda):
opt_t_1 := Optimization:-NLPSolve(op(1, v), t=0..Pi, maximize):
opt_t_2 := Optimization:-NLPSolve(op(1, v), t=Pi..2*Pi, maximize):

opt_x_1 := op(2, v) *~ eval([cos(t), sin(t)], opt_t_1[2]):
opt_x_2 := op(2, v) *~ eval([cos(t), sin(t)], opt_t_2[2]):

identify(evalf[8](opt_x_1));
identify(evalf[8](opt_x_2));

              [1  (1/2)       2  1  (1/2)       2]
              [- 2      lambda , - 2      lambda ]
              [2                 2               ]

            [  1  (1/2)       2    1  (1/2)       2]
            [- - 2      lambda , - - 2      lambda ]
            [  2                   2               ]

Then the maxima are rejected to x = +/- infinity and y = x.

Which s a far better result than Mathematica provides.

So, let me ask you this question (and @QM too): suppose your are part of a MMA users blog and MMA tells you that

 {\[Infinity], {x -> Indeterminate, y -> Indeterminate}}

wouldn't you say that MMA doesn't give a completely satisfactory answer ?

I feel a sense of unfairness in the question and this two-man bashing of Maple, am I wrong?

d1 is a centered divided difference appr...

Answered: mmcdara 7891

March 20 2023

1 2

d1 is a centered divided difference approximation of diff(u, x$3) at "point" [i+1/2, j];not an approximation of the laplacian of u at [i, j]

>

restart

>	tau := (d, s) -> p -> p + s*~[2-d, d-1]

(1)

>	# examples tau(1, h)([0, 0]), tau(1, -h)([0, 0]); tau(2, k)([0, 0]), tau(2, -k)([0, 0]);

(2)

>	# forward and backward divided differences as approximations of diff(u(x, y), x) fd_x := (u[op(tau(1, 1)([i, j]))] - u[op(tau(1, 0)([i, j]))]) / h; bd_x := (u[op(tau(1, 0)([i, j]))] - u[op(tau(1, -1)([i, j]))]) / h;

(3)

>	# approximation of diff(u(x, y), x$2) from fd and bd d2_x := simplify((fd_x-bd_x)/h)

-(-u[1+i, j]+2*u[i, j]-u[-1+i, j])/h^2

(4)

>	# laplacian at point [i, j] fd_y := (u[op(tau(2, 1)([i, j]))] - u[op(tau(2, 0)([i, j]))]) / k; bd_y := (u[op(tau(2, 0)([i, j]))] - u[op(tau(2, -1)([i, j]))]) / k; d2_y := simplify((fd_y-bd_y)/k); lap := simplify(d2_x+d2_y): lap := collect(numer(lap), u[i, j]) / denom(lap)

-(-u[i, 1+j]+2*u[i, j]-u[i, -1+j])/k^2

((-2*h^2-2*k^2)*u[i, j]+h^2*u[i, -1+j]+h^2*u[i, 1+j]+k^2*u[-1+i, j]+k^2*u[1+i, j])/(h^2*k^2)

(5)

>	# approximation of diff(u(x, y), x$3) at point [i+1/2, j] d3_x := simplify((eval(d2_x, i=i+1) - d2_x) / h)

(u[2+i, j]-3*u[1+i, j]+3*u[i, j]-u[-1+i, j])/h^3

(6)

>

Download DividedDifferences.mw

This code enables printing your list in ...

Answered: mmcdara 7891

March 19 2023

1 3

This code allows you to print your list in a way that is close to what you want.
Printed lines have arround n=10 characters.

Of course arrangements have to be done if some elements of the list are numbers with more than n digits.
More of that, the list is converted into a string before being parsed, which means that spaces and commas are included in the calculation of the string to plot.
But I think this can give you some hints to go further (for instance if you want to plot numbers, you must augment the the range to plot 1..pc[cut] by 2*p, where p is number of occurrences of a comma in newLs[1..pc[cut]].

>

restart:

>

SmartPrint := proc(L::list, n::posint)
local newLs, pc, cut:

newLs := convert(L, string)[2..-2]:
pc     := [ StringTools:-SearchAll(",", newLs)]:
cut    := ListTools:-BinaryPlace(pc, n+1/2):

while length(newLs) > n do
  if pc[cut+1] = n+1 then
    cut := cut+1:
  end if:
  printf("%-20s (length = %d)\n", StringTools:-Trim(newLs[1..pc[cut]]), pc[cut]):
  newLs := StringTools:-Trim(newLs[pc[cut]+1..-1]);
  pc     := [ StringTools:-SearchAll(",", newLs)];
  cut    := ListTools:-BinaryPlace(pc, n+1/2);
end do:

printf("%-20s (length = %d)\n", StringTools:-Trim(newLs), length(StringTools:-Trim(newLs))):
end proc:

>	L := [seq(i,i=1..20)]: SmartPrint(L, 10)

1, 2, 3, 4,          (length = 11)
5, 6, 7, 8,          (length = 11)
9, 10, 11,           (length = 10)
12, 13, 14,          (length = 11)
15, 16, 17,          (length = 11)
18, 19, 20           (length = 10)

>	r := 1000rand(0..1)rand(0..1) + 100rand(0..1)rand(0..1) + 10rand(0..10)rand(0..1) + rand(0..10): L := [seq(r(), i=1..20)]; SmartPrint(L, 10)

60, 1100,            (length = 9)
13, 1, 110,          (length = 11)
66, 2, 2,            (length = 9)
101, 1008,           (length = 10)
28, 1038,            (length = 9)
100, 3, 8,           (length = 10)
33, 1161,            (length = 9)
4, 50, 101           (length = 10)

>

Download SmartPrint.mw

Either you miss one ODE or you need to p...

Answered: mmcdara 7891

March 18 2023

1 1

Unless I'm mistaken you miss one equation:

eq1 is of the form
```
 A*diff(h(r), r$4) + B*diff(h(r), r$4) + LowerOderTerms = 0
```
A and B do not contain derivatives of h anq q of order 4
eq2 is of the form
```
 C*diff(h(r), r$3) + E*diff(h(r), r$3) + LowerOderTerms = 0
```
C and E do not contain derivatives of h anq q of order larger than 2
Then plugging eq2 into eq1 produces a single ODE of the formeq1 is of the form
```
 P*diff(h(r), r$4) + Q*diff(h(r), r$4) + LowerOderTerms = 0
```
where P and Q do not contain derivatives of h anq q of order larger than 3

So you miss one ODE.
The detailed explanation is given in the attached file.
A toy problem is also provided to help you understand why dsolve doesn't work asyou expected.

>

restart

>

(1)

>

Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

>	# check if the system {eq1, eq2} is an ODE system wrt h(t) and q(t); sol := solve({eq1, eq2}, [diff(h(r), r$4), diff(q(r), r$4)])

(2)

>	# as it's not the case try to represent rhis system in a simple form with(LargeExpressions):

>	H1 := collect(eq1, diff(h(r), r$4), Veil[A] ); # check if A[1] contains diff(q(r), r$4) Unveil[A](A[1]); # Veil A[2] Q1 := collect(expand(Unveil[A](A[2])), diff(q(r), r$4), Veil[A2] );

(2/45)*A[1]*(diff(diff(diff(diff(h(r), r), r), r), r))+(1/360)*A[2] = 0

(r*(diff(q(r), r))+2*q(r))*q(r)^2*(r^2*q(r)+h(r))^2*h(r)^2/(Pi*(-r*q(r)*(diff(h(r), r))+2*q(r)*h(r)+r*(diff(q(r), r))*h(r))^5)

16*A2[1]*(diff(diff(diff(diff(q(r), r), r), r), r))+A2[2]

(3)

>	# Then eq1 is of the form edo1 := ''eq1'' = eval(H1, A[2]=Q1): edo1;

eq1 = ((2/45)*A[1]*(diff(diff(diff(diff(h(r), r), r), r), r))+(2/45)*A2[1]*(diff(diff(diff(diff(q(r), r), r), r), r))+(1/360)*A2[2] = 0)

(4)

>	# observe that eq2 doen't contain any derivatives of order 4 indets(eq2);

{r, diff(diff(diff(h(r), r), r), r), diff(diff(diff(q(r), r), r), r), diff(diff(h(r), r), r), diff(diff(q(r), r), r), diff(h(r), r), diff(q(r), r), h(r), q(r)}

(5)

>	H2 := collect(eq2, diff(h(r), r$3), Veil[B] ); # check if B[1] contains diff(q(r), r$3) Unveil[B](B[1]); # Veil B[2] Q2 := collect(expand(Unveil[B](B[2])), diff(q(r), r$3), Veil[B2] ); # Then eq2 is of the form edo2 := ''eq2'' = eval(H2, B[2]=Q2): edo2;

(1/45)*B[1]*(diff(diff(diff(h(r), r), r), r))+(1/360)*B[2] = 0

h(r)*(r^2*q(r)+h(r))^2*(2*r*(diff(q(r), r))*h(r)+r*q(r)*(diff(h(r), r))+4*q(r)*h(r))*q(r)*(r*(diff(q(r), r))+2*q(r))/((-r*q(r)*(diff(h(r), r))+2*q(r)*h(r)+r*(diff(q(r), r))*h(r))^5*r*Pi)

8*B2[1]*(diff(diff(diff(q(r), r), r), r))+B2[2]

eq2 = ((1/45)*B[1]*(diff(diff(diff(h(r), r), r), r))+(1/45)*B2[1]*(diff(diff(diff(q(r), r), r), r))+(1/360)*B2[2] = 0)

(6)

>	# thus your system writes: # (A, A2, are functions of r and derivatives of f and q up to order 3, at most) # (B, B2, are functions of r and derivatives of f and q up to order 2, at most) print~({edo1, edo2}):

eq1 = ((2/45)*A[1]*(diff(diff(diff(diff(h(r), r), r), r), r))+(2/45)*A2[1]*(diff(diff(diff(diff(q(r), r), r), r), r))+(1/360)*A2[2] = 0)

eq2 = ((1/45)*B[1]*(diff(diff(diff(h(r), r), r), r))+(1/45)*B2[1]*(diff(diff(diff(q(r), r), r), r))+(1/360)*B2[2] = 0)

(7)

>	# Veil A2 wrt diff(h(r), r$3) collect(expand(Unveil[A2](A2[2])), diff(h(r), r$3), Veil[A22] ): collect(expand(Unveil[A2](A2[3])), diff(h(r), r$3), Veil[A23] ): # thus eq1 writes edo1_1 := ''eq1'' = eval(rhs(edo1), [A2[2]=%%, A2[3]=%]): edo1_1;

eq1 = ((2/45)*A[1]*(diff(diff(diff(diff(h(r), r), r), r), r))+(2/45)*A2[1]*(diff(diff(diff(diff(q(r), r), r), r), r))+(1/45)*A22[1]*(diff(diff(diff(h(r), r), r), r))+(1/360)*A22[2] = 0)

(8)

>

# If you plug the expression of eq2 into the expression above, you finally get a single ODE
# with two unknowns and then you miss one ODE.
#
#
# To convince you look at thimple toy model

edosys := {
  diff(f(x), x$2)+diff(g(x), x$2)=x,
  diff(f(x), x)+diff(g(x), x)=1,
  f(0)=0, D(f)(0)=0,
  g(0)=0, D(g)(0)=0
};
infolevel[dsolve]:=100:
dsolve(edosys);

{diff(diff(f(x), x), x)+diff(diff(g(x), x), x) = x, diff(f(x), x)+diff(g(x), x) = 1, f(0) = 0, g(0) = 0, (D(f))(0) = 0, (D(g))(0) = 0}

(9)

>	# and now to this well constructed toy model edosys := { diff(f(x), x$2)+diff(g(x), x$2)=x, diff(f(x), x$2)-diff(g(x), x$2)=1, f(0)=0, D(f)(0)=0, g(0)=0, D(g)(0)=0 }; dsolve(edosys);

{diff(diff(f(x), x), x)-(diff(diff(g(x), x), x)) = 1, diff(diff(f(x), x), x)+diff(diff(g(x), x), x) = x, f(0) = 0, g(0) = 0, (D(f))(0) = 0, (D(g))(0) = 0}

-> Solving each unknown as a function of the next ones using the order: [g(x), f(x)]
-> Calling odsolve with the ODE diff(diff(y(x) x) x) = -1/2+(1/2)*x y(x) explicit singsol = none
Methods for second order ODEs:
--- Trying classification methods ---
trying a quadrature
<- quadrature successful
-> Calling odsolve with the ODE diff(diff(y(x) x) x) = 1/2+(1/2)*x y(x) explicit singsol = none
Methods for second order ODEs:
--- Trying classification methods ---
trying a quadrature
<- quadrature successful

${f(x) = (1/12)*x^3+(1/4)*x^2, g(x) = (1/12)*x^3-(1/4)*x^2}$

(10)

>

Download dsolve_mmcdara.mw

A toy_problem dsolve can't find the solution of, but whose solution can be obtained after reparameterisation:

>

restart:

>	indolevel[dsolve] := 100:

>	# A toy problem edosys := { diff(f(x), x$2)+diff(g(x), x$2)=0, diff(f(x), x)+diff(g(x), x)=0, f(0)=0, D(f)(0)=1, g(0)=0, D(g)(0)=-1 }; infolevel[dsolve]:=100: sol := dsolve(edosys);

{diff(diff(f(x), x), x)+diff(diff(g(x), x), x) = 0, diff(f(x), x)+diff(g(x), x) = 0, f(0) = 0, g(0) = 0, (D(f))(0) = 1, (D(g))(0) = -1}

-> Solving each unknown as a function of the next ones using the order: [g(x), f(x)]
-> Calling odsolve with the ODE diff(y(x) x) = -(diff(f(x) x)) y(x) explicit singsol = none

Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature

<- quadrature successful

>	# This system writes edosys := [ diff(df(x), x)+diff(dg(x), x)=0, diff(f(x), x)=df(x), diff(g(x), x)=dg(x), df(x)+dg(x)=0, f(0)=0, df(0)=1, g(0)=0, dg(0)=-1 ]: print~(edosys): sol := dsolve(edosys);

-> Solving each unknown as a function of the next ones using the order: [g(x), f(x), dg(x), df(x)]
-> Calling odsolve with the ODE diff(y(x) x) = -df(x) y(x) explicit singsol = none
Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature

<- quadrature successful
-> Calling odsolve with the ODE diff(y(x) x) = df(x) y(x) explicit singsol = none
Methods for first order ODEs:
--- Trying classification methods ---
trying a quadrature
<- quadrature successful

>	# Observe that original unknowns f(x) and g(x) are meaningless edosys := [ diff(df(x), x)+diff(dg(x), x)=0, df(x)+dg(x)=0, df(0)=1, dg(0)=-1 ]: print~(edosys): sol := dsolve(edosys);

(1)

>

Download Toy_Problem.mw

infolevel[solve] may help...

Answered: mmcdara 7891

March 18 2023

1 3

restart:
sys := {x^2*y+x=1, x*y-y^2=0}:
solve(sys);
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  / 
infolevel[solve]:=100:
solve(sys)
Main: Entering solver with 2 equations in 2 variables
Main: attempting to solve as a linear system
Main: attempting to solve as a polynomial system
Main: Polynomial solver successful. Exiting solver returning 1 solution
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  /

This suggests solve calls some polynomial solver.
Searching for it in the help pages returns SolveTools/PolynomialSystem. I don't know if solve calls this procedure (showstat(solve) is hard to browse). Assuming it does:

SolveTools:-PolynomialSystem(sys, {x,y});

Main: polynomial system split into 1 parts under preprocessing
Main: trying resultant methods
PseudoResultant: normalize equations, length= 43
PseudoResultant: 353502 [1 3999999986 y 1] 2 2 43 0 3 0
PseudoResultant: normalize equations, length= 13
PseudoResultant: normalize equations, length= 35
PseudoResultant: 30000 [1 100000559 x] 1 1 13 0 3 1
PseudoResultant: normalize equations, length= 3
PseudoResultant: -10 [] 0 0 3 0 3 1
PseudoResultant: 201729 [1 1000012108 y] 2 2 35 1 8 0
PseudoResultant: normalize equations, length= 21
PseudoResultant: 117197 [1 700003740 x] 1 1 21 1 8 0
PseudoResultant: normalize equations, length= 3
PseudoResultant: -10 [] 0 0 3 0 3 0
PseudoResultant: 2 solutions found, now doing backsubstitution
PseudoResultant: backsubstitution of x
PseudoResultant: backsubstitution of y
PseudoResultant: backsubstitution of x
PseudoResultant: backsubstitution of y
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  /

The simplest answer is...

Answered: mmcdara 7891

March 17 2023

1 0

restart
p := 35*x^8 + 80*x^7 + 140*x^6 + 504*x^5 + 1260*x^4 + 1960*x^3 + 1960*x^2 + 1120*x + 10:
dp := diff(p, x):
fsolve(dp);
                         -1., -1., -1.

But I suspect that you are looking for an answer you can obtain "by hand" aren't you?
Here is, written in Maple, a step-by-step procedure which only requires doing polynomial division (which can be done by hand)
Download ByHand.mw

Details:

Observe -1 is a root of dp = diff(p, x)
Thus dp = Pol(x) * (x+1)^r with r >=1 and degree(Pol) <= 6
Find the value of r this way:
1. Set Q =Pol(x)
2. Divide (can be done by hand) Q by (x+1) and note R the remainder
3. Go to point 1 until R = 0
The value of r is the number of times you executed the steps 1 to 3
Here r=3 ==> dp = Pol(x) * (x+1)^3 with degree(Pol) = 4
Search if Pol(x) has a real zero Z by applying the procedure given at points 1 to 3.
Here you recursively divide 4 times Pol(x) by (x-Z) : if the last remainder is 0, then is a root of dp, if not (which is the case) dp has only 1 real root (-1) of multipliciy 3.

Thus p has a only one stationnary point at x=-1 (I don't know if you must answer 1 or 3 ???)

One way...

Answered: mmcdara 7891

March 17 2023

1 1

The goal is to produce almost exactly the expressio of the solution you give, thus we have to go beyond the dsolve to reshape the solution Maple provides.
The procedure is described step-by-step

>

restart:

>

local I: