sand15

The smallest value of Em is......

Answered: sand15 1182

Yesterday at 4:20 PM

1 5

-0.06888275721

You write "I want to choose A to obtain the smallest absolute value of Em, preferably ~ -0.00005".
I suppose you want to say "I want to choose A to obtain the smallest absolute value of Em, preferably ~ +0.00005", don't you?

Whatever, Em is strctly negative for A > 0 and its limit when A goes to 0 from the right is -0.06888275721, which means that your requirement cannot be met.

Note that for A = 0 we have E(𝜐) = C𝜐² + -0.06888275721 which immediatly gives the value of C... but the fitted model is unacceptable has shown in the attached file.

At the very end I'm leaning to consider your problem is not correctly posed and I suggest you to consider carefully what you want to do.
ML_sand15.mw

Huge errors...

Answered: sand15 1182

October 31 2025

1 7

for i2 varying from 0 to 1

doen't nmean anything at all as the
Error, controlling variable of for loop must be a name or sequence of 2 names
says: between for and from you must have only one name (here you have i2 and varying).

Given the past worksheets you posted I don't even understand how you can have written (sorry if I feel that rude) such a stupid thing?
A few correct syntaxes are

for i2 from 0 to 1 do
...
end do:

for i2 from 0 to 1 by 1/3 do
...
end do:

# etc

Note that the first syntax seems quite stupid to me in your context:

# with K2 = 0.1968052257

for i2 from 0 to 1 do
  if i2 > K2 then 
    do_something
  else
    do_something_else
  endif:
end do:

because the counter (here i2) being an integer unless "by some_value" is specified, it takes only the values i1=0 and i2=1. So a simpler operation is

Do_this := i2 -> piecewise(i2 > K2, do_something, do_something_else)

Nothing of what you try to achieve is clear to me.
Moreover, if you are not capable to manage correctly the 2D input document mode, just as I can't neither, use the 1D input worksheet mode (your file contains a lot of breaks which generate errors (yout for anf if structures do not end correctly).

Whatever... look to the attached file to see how to write something error free A_no_error_worksheet_sand15.mw and try to use it contents to formulate correctly your problem (which I repeat I do not understand).

One way...

Answered: sand15 1182

October 31 2025

3 0

(done with Maple 2015 but should work the same with Maple 2018)
One_way_sand15.mw (read carefully the last comment in thiw worksheet)

Proceed the same way for the 3D plots by replacing p1_final by p1 in C11 and p2_final by p2 in C12.

Note: If you do some browsing on this same site you should be able to find several questions about the same subject.
@acer already gave a lot of answers too [moderator: eg. here, here, here]
[moderator: If you need to export the 3D plot with colorbar in such an old Maple version then see also here.]

[moderator: colorbars received direct support for 3D plotting commands in Maple 2024.]

One way...

Answered: sand15 1182

October 30 2025

2 1

restart;
n := 0:
L := []:

for a from 3 to 100 do
    for b from 3 to a do
        c2 := a^2 - a*b + b^2;
        c := isqrt(c2);
        if c^2 = c2 then
            if c < a + b and a < b + c and b < c + a then
                if igcd(a, b, c) = 1 then
                    x2 := (-a^2 + b^2 + c^2)/2;
                    y2 := (a^2 - b^2 + c^2)/2;
                    z2 := (a^2 + b^2 - c^2)/2;
                    if 0 < x2 and 0 < y2 and 0 < z2 then
                        # One way
                        S := [sqrt(x2), sqrt(y2), sqrt(z2)];
                        S := S[sort(evalf(S), output=permutation)];
                        x := S[1]; 
                        y := S[2]; 
                        z := S[3];
                        n := n + 1;
                        L := [op(L), [x, y, z, sqrt(x^2 + y^2), sqrt(y^2 + z^2), sqrt(x^2 + z^2)]];
                    end if;
                end if;
            end if;
        end if;
    end do;
end do;

n;
L;

Why was sort([sqrt(x2), sqrt(y2), sqrt(z2)]) ineffective?
Look at this

S0 := [6*sqrt(110), 2*sqrt(610), 3*sqrt(649)];
sort(S0, `>`)
Error, invalid input: a numeric list is expected for sort method ``>``

# To understand this "failure" look at this 
if S0[1] > S0[2] then "yes" else "no" end if;
Error, cannot determine if this expression is true or false: 2*610^(1/2) < 6*110^(1/2)

So sort, which vasically uses an if structure, simply doesn't know how to compare these two non rational numbers.

Note that is can do it:

is(S0[1] > S0[2]);
                             true
is(S0[1] < S0[2]);
                             false

Another way...

Answered: sand15 1182

October 26 2025

2 5

(strangely it takes a long time to be processed)

>

restart

>	kernelopts(version)

(1)

>	expr := sqrt((-alpha^2x^2 + R^2 + x^2)/(-alpha^2x^2 + R^2))/sqrt(R^4/((-alpha^2x^2 + R^2 + x^2)^2(-alpha^2*x^2 + R^2))):

>	subexpr := indets(expr, `+`); CodeTools:-Usage( simplify(expr) assuming op(subexpr >=~ 0), R >= 0 )

${-alpha^2*x^2+R^2, -alpha^2*x^2+R^2+x^2}$

memory used=0.59GiB, alloc change=429.34MiB, cpu time=8.33s, real time=8.16s, gc time=356.09ms

(-alpha^2*x^2+R^2+x^2)^(3/2)/R^2

(2)

Download simplify_radical_sand15.mw

One way...

Answered: sand15 1182

October 17 2025

1 1

From the code below you should be able to define any binary operator you want.

>

restart

>	kernelopts(version)

(1)

>	Isom := `&cong;`: Isom(A, B);

(2)

>	IsContained := `⊆`: IsContained(A, B)

(3)

>	IsNotContained := `⊈`: IsNotContained(A, B)

(4)

>	DoesNotContain := `⊉`: DoesNotContain(A, B)

(5)

>	# For details look to the answers to this QUESTION `print/Fs` := proc(x,y) uses Typesetting; mrow( Typeset(x), mo("⁢"), mover(mo("→"),mo("s")), mo("⁢"), Typeset(y) ); end proc: Fs(A, B)

(6)

>	# For any other hexadecimal code see HERE , for instance here are different arrows: `#mo("&mapsto;")`; `#mo("↣")`; `#mo("↠")`;

(7)

Download Special_Operators.mw

value...

Answered: sand15 1182

October 12 2025

1 2

AA := value( eval(A, {Omega = (t -> 1), chi = (t -> 1), g = (t -> 1), i = 1, p = (t -> 1), i = 1}) );
has(AA, {int, Int})
                             false

h-pde-M_sand15.mw

Something like that?...

Answered: sand15 1182

October 07 2025

1 1

t-plot_sand15.mw

Some functions, like 'sin' for instance, leads to an unevaluated integral.

Even using evalf/Int doesn't help for the the numeric values are Float(undefined).

Something like this?...

Answered: sand15 1182

October 07 2025

2 7

Download DGL_test_sand15.mw

Explanations...

Answered: sand15 1182

October 03 2025

2 6

@JoyDivisionMan

Using "blindly" Maple (2015) provides a wrong result (see attached file) for, I guess, maple does not handle correctly the absolute value.

So here is he trick.
Let f_Z(z) the pdf of Z = |Y₁- Y₂| and f_U(u) the pdf of U = Y₁- Y₂.
Of course the support of Z is (0, 2) while U's is (-2, 2). But Y₁ and Y₂ both having the same symmetic distribution it is easy to see that f_U(u) is symmetric too.
So consider only the restriction of f_U(u) the its half support (0, 2) and twice this latter: then 2f_U(u) = f_Z(u).

With that trixk Maple provides the corret result.

Note that f_U(u) can be computed using a convolution product.
Indeed U = Y₁ - Y₂ = Y₁ + (-Y₂) = Y₁ + Y₂ because f_Y(y) is symmetrc wrt y=0.
Then f_U(u) = (f_Y * f_Y)(y) where '*' denotes the convolution produce.

Please, if that's okay with you, don't forget (as it was the case with your last question) to let me know.
(Doing work for someone without knowing whether they appreciated it or not is never pleasant.)

>

restart

>	kernelopts(version)

(1)

>	with(Statistics):

>	f := x -> piecewise(abs(x) < 1, 2*sqrt(1 - x^2)/Pi, 0);

proc (x) options operator, arrow; piecewise(abs(x) < 1, 2*sqrt(1-x^2)/Pi, 0) end proc

(2)

>	F := unapply(int(f(t), t=-1..x), x);

proc (x) options operator, arrow; piecewise(x <= -1, 0, x <= 1, (1/2)*(2*x*(1-x^2)^(1/2)+2*arcsin(x)+Pi)/Pi, 1 < x, 1) end proc

(3)

>	d := Distribution( PDF = (x -> f(x)), CDF = (x -> F(x)) ):

>	Y__1 := RandomVariable(d): Y__2 := RandomVariable(d):

>	# From symmetry considerations the PDF of Z is twice the PDF of U restricted to U > 0 Z := abs(U): U := Y__1 - Y__2:

>	# Note the ∞ term in the case 0 < z < 1 !!! # This is obviously a bug. PDF(Z, z);

(4)

>	# To help Maple find the correct PDF we may observe that, from symmetry considerations, # the PDF of Z is twice the PDF of U = Y__1 - Y__2 once restricted to U > 0 U := Y__1 - Y__2: f__U := unapply( (2*PDF(U, u) assuming u > 0), u);

(5)

>	S__Z := Sample(Z, 10^6): plots:-display( Histogram(S__Z, transparency=0.7, minbins=100, style=polygon) , plot(f__U(u), u=0..2, color=red, thickness=2, legend=typeset(:-PDF('Z'))) )

>

F__U := proc(u)
  option remember:
  local F1;
  F1 := 2*evalf(Int(op([2, 2], f__U(t)), t=0..1)):
  if _passed::numeric then
    if u <= 1 then
      2*evalf(Int(op([2, 2], f__U(t)), t=0..u))
    elif u <= 2 then
      F1 + 2*evalf(Int(Re(op([2, 4], f__U(t))), t=1..u, method = _d01ajc))
    end if:
  else
    procname(_passed)
  end if:
end proc:

>	plot([seq([u, F__U(u)], u=0..2, 0.1)], color=red, thickness=2, legend=typeset(:-CDF('Z')), gridlines=true)

>	# Because Y__1 and Y__2 both have the same symmetric PDF: # # U = Y__1 - Y__2 = Y__1 + (-Y__2) = Y__1 + Y__2 # # Thus the PDF of U can be obtained through a convolution product f__Y := unapply(PDF(Y__1, y), y): C := unapply( int(f__Y(y)*f__Y(u-y), y = -infinity..+infinity), u ):

>	plots:-display( Histogram(S__Z, transparency=0.7, minbins=100, style=polygon) , plot(2*C(u), u=0..2, color=red, thickness=2, legend=typeset(:-PDF('Z'))) )

Depending on what statistic you want it is sometimes more powerfull to use U instead of Z

>	:-Mean('Z') = Mean(Z);

Mean(Z) = (256/45)/Pi^2

(6)

>	:-Variance('Z') = Variance(Z);

Variance(Z) = (1/4050)*(2025*Pi^4-131072)/Pi^4

(7)

>	:-Median('Z') = Median(Z); :-Median('Z') = Median(Z, numeric);

(8)

>	:-Median('Z') = fsolve('F__U'(u)=0.5, u=0..2); med := rhs(%):

(9)

Spline approximation of F__U

>	S__U := unapply(CurveFitting:-Spline([seq([u, F__U(u)], u=0..2, 0.25)], u), u):

Procedure to compute quantiles of Z

>

quantile := proc(p)
  local guess:
  if F__U(1) > p then
  # guess := fsolve(S__U(u)=p, u=0..1)
    guess := fsolve('F__U'(u)=p, u=0..1)
  else
    guess := fsolve(S__U(u)=p, u=1..2)
  end if:
  return [guess, ``(F__U(guess))]:
end proc:

:-Quantile('Z', 0.05) = quantile(0.05); # Exact, based on F__U
:-Quantile('Z', 0.80) = quantile(0.80); # Exact, based on F__U
:-Quantile('Z', 0.85) = quantile(0.85); # Approximation based on S__U
:-Quantile('Z', 0.95) = quantile(0.95); # Approximation bBased on S__U

(10)

>	L := limit(f__U(u), u=0, right); if false then plots:-display( plot(f__U(u), u=0..2, color=red), plottools:-transform((x, y) -> [2-x, L-y])(plot(f__U(u), u=0..2, color=blue)) ); end if;

(32/3)/Pi^2

(11)

Download PDF_Z_sand15.mw

KKT doesn't seem necessary here...

Answered: sand15 1182

September 30 2025

1 12

𝛑_m is a priori quadratic wrt i₁ so 𝛑_m has only one extremum swhose status (minimum vs maximum) depends on the sign of the coefficient of i₁².

Note: I say "𝛑_m is a priori quadratic wrt i₁" for 𝛑_m writes P . i₁² + Q . i₁ + R where P, Q, R are expressions which depend on several parameters. It could therefore happen that, depending on the values of these parameters, P, or even P and Q, might be equal to 0.
The case P = 0, Q <> 0 is even simpler because 𝛑_m necessary reaches its the minimum at one of the two boundaries defined by A and B (see below the meaning of these two quantities)

As you did not provide more informations I considered arbitrarily that P was different from 0.

Let x* the location of this extremum and D the domain defined by the two primal feasibility conditions

f₁ : A <= i₁
f₂ : i₁ <= B (B supposedly larger than A of course)

Assuming for instance that P <> 0 and that 𝛑_m reaches its minimum at some point i₁= x*, three cases are to be examined:

x* < A:
The constrained minimum is x* = A
x* > B:
The constrained minimum is x* = B
A < x* > B:
The constrained minimum is x*

Assuming now that P = 0 and Q <> 0, we find the location of the minimizer i₁= x* of 𝛑_m corresponds to one of these two cases:

𝛑_m(A) < 𝛑_m(B)
The constrained minimum is x* = A
𝛑_m(A) > 𝛑_m(B)
The constrained minimum is x* = B

I don't really think it's necessary to write KKT conditions and solve grad(Lagrangian(i₁, 𝛍₁, 𝛍₂)) = 0 wrt i₁, 𝛍₁, 𝛍₂ to get this result.

Basically your problem can be parameterized by the four quantities P, Q, A, B.

Why not simply use solve/parametric to get all the possible solutions?

>

restart

>	SC_CS_sols := solve( { diff(Pi1^2 + Qi1 + R, i1) - mu[1] + mu[2] = 0, mu[1](A-i1) = 0, mu[2](i1-B) = 0 } , {i1, mu[1], mu[2]} , 'parametric'='full', 'parameters'={P, Q, A, B} )

-

(1)

>	pi__m := (w-i1)(1/2+(i1-i2)/(2tau))(1-(Pn-Pr)/(1-delta))-C0-(1/2)Cr(1/2+(i1-i2)/(2tau))^2(1-(Pn-Pr)/(1-delta))^2+Cerho0*(1-(Pn-Pr)/(1-delta)): collect(pi__m, i1): PQR := [R, Q, P] =~ [coeffs(%, i1)]: print~(PQR):

R = w*(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))-C0-(1/2)*Cr*(1/2-(1/2)*i2/tau)^2*(1-(Pn-Pr)/(1-delta))^2+Ce*rho0*(1-(Pn-Pr)/(1-delta))

Q = ((1/2)*w/tau-1/2+(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))-(1/2)*Cr*(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))^2/tau

P = -(1/2)*(1-(Pn-Pr)/(1-delta))/tau-(1/8)*Cr*(1-(Pn-Pr)/(1-delta))^2/tau^2

(2)

>	AB := [A = (2rho0-1)tau+i2, B = tau+i2]

(3)

>	ind := indets(rhs~(PQR)) union indets(rhs~(AB));

(4)

>

randomize();

data := { seq(i = rand(0. .. 1.)(), i in ind) };
PQRAB := { eval(PQR, data)[], eval(AB, data)[] };

all_solutions := allvalues~(eval(SC_CS_sols, PQRAB)):
all_pi        := map(s -> eval(eval(pi__m, data), s), all_solutions):
here          := ListTools:-Search(min(all_pi), all_pi):
SOLUTION      := [all_pi[here], all_solutions[here]];

plot(eval(pi__m, data), i1=eval(eval(A..B, AB), data), gridlines=true)

>

Download proposal_sand15.mw

elementwise comparisons...

Answered: sand15 1182

September 28 2025

0 2

See the Maple help page for the meaning of the elementwise operator '~'

a := Matrix([[1, 2], [1, 2]]):
is~(a[1] =~ a[2]);
    
b := [seq(a[i], i=1..2)]:
is~(b[1] =~ b[2]);

Suggestion...

Answered: sand15 1182

September 28 2025

0 1

Look HERE to my reply to a recent similar question.

Simply...

Answered: sand15 1182

September 27 2025

3 2

The lazzy attitude: simply ask Maple the result

>

restart

>	with(Statistics):

>	f[1] := t -> piecewise(t <= 0, 0, 0 < t and t < 2, 1/(Pisqrt(1-(1-t)^2)), t >= 2, 0); f[2] := t -> piecewise(t <= 0, 0, 0 < t and t < 1, 2t, t >= 1, 0);

proc (t) options operator, arrow; piecewise(t <= 0, 0, 0 < t and t < 2, 1/(Pi*sqrt(1-(1-t)^2)), 2 <= t, 0) end proc

(1)

>	Dist1 := Distribution(PDF = (t -> f[1](t))): Dist2 := Distribution(PDF = (t -> f[2](t))):

>	X := RandomVariable(Dist1): Y := RandomVariable(Dist2): Z := X*Y:

>	`PDF(Z)` := PDF(Z, t); `CDF(Z)` := CDF(Z, t);

`PDF(Z)` := piecewise(t < 0, 0, t <= 2, -(2/3)*(t^2-t-2)/(sqrt(t)*sqrt(2-t)*Pi), 2 < t, 0)

piecewise(t <= 0, 0, t < 2, (1/6)*(2*t^(3/2)*(2-t)^(1/2)*(-t*(t-2))^(1/2)+2*t^(1/2)*(2-t)^(1/2)*(-t*(t-2))^(1/2)+6*t^(1/2)*(2-t)^(1/2)*arcsin(-1+t)+3*Pi*(-t*(t-2))^(1/2))/(Pi*(-t*(t-2))^(1/2)), 2 <= t, 1)

(2)

>	plot(`PDF(Z)`, t=0..2)

>

Download Inquiry_sand15.mw

A classical result (under some restricions)
Let X and Y two independent random variables whose both supports belong to (0, +∞) and with respective pdf f_X and f_Y then the pdf f_Z of Z = X Y is given by

For more details you can see for instance Glen at al.or go to Google Scholar and search for "product of two independent random variables".

If you want to understand how to get the result while not knowing about this formula...
As a rule
To obtain the probability density function of the function h of one (or more) continuous random variable(s) it is always simpler to compute its cumulative function first and next to differentiate it (assuming this cumulative function is derivable).
Thus is all the more simple if h has an unique inverse over the support of this (these) random variable(s) (if not some precautions have to be taken).
This your case: Z = h( X , Y ) = X Y.
Worksheet PROCEDURE.mw describes the step-by-step procedure to get f_Z :

>

restart

Let X and Y two independent random variables with respective probability functions

>	f__X := x -> piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pisqrt(1-(1-x)^2)), x >= 2, 0); f__Y := y -> piecewise(y <= 0, 0, 0 < t and y < 1, 2y, y >= 1, 0);

proc (x) options operator, arrow; piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0) end proc

(1)

Their cumulative functions write

>	F__X := unapply(int(f__X(t), t=-infinity..x), x); F__Y := unapply(int(f__Y(t), t=-infinity..y), y);

(2)

Let C = X Y.

>

Prob(Z < z) = Prob(XY < z);
# Independency of X and Y gives

Prob(Z < z) = Prob(X < x and Y < z/x);

# Prob(Z < z | X < x) reads: Prob(Z < z) given that X < x

Prob(cat(Z < z, ` | `, X = x)) = Prob(Y < z/x);

# By definition of F_Y

Prob(cat(Z < z, ` | `, X = x)) = F__Y(z/x);

# To get Prob(Z < z) integrate Prob(Z < z | X < x) for all x values (accounting for the PDF of X of course)
Prob(Z < z) = Int(F__Y(z/x)*f__X(x), x=0..+infinity);

# By definition of F_Z

F__Z(z) = rhs(%);

# Derivation both sides wrt z gives the formal expression of f_Z

f__Z(z) = Diff(rhs(%), z);

# Exchange Diff and Int (assuming this operation is licit)

f__Z(z) = Int(Diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# Evaluate the derivative of the integrand

f__Z(z) = Int(diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# Evaluate the integral of the previous result

f__Z(z) = int(diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# (You see that the Mellin transform is absolutely useless here.)

Prob(`Z < z | X = x`) = piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)

Prob(Z < z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

`#msub(mi("F"),mi("Z"))`(z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

`#msub(mi("f"),mi("Z"))`(z) = Diff(Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity), z)

`#msub(mi("f"),mi("Z"))`(z) = Int(Diff(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), z), x = 0 .. infinity)

`#msub(mi("f"),mi("Z"))`(z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, 2*z/x^2, 1 < z/x, 0)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

f__Z(z) = piecewise(z < 0, 0, z < 2, -(2/3)*(z^2-z-2)/((z*(2-z))^(1/2)*Pi), 2 <= z, 0)

(3)

>	plot(rhs(%), z=0..2)

Worksheet PROCEDURE_formal.mw contains the demonstration of the classical result when X and Y are independent random variables with supports in (0, +∞).

As you already asked a question about the convolution product, here is a useful trick
Let A=ln( X ) and B=ln( Y ), the pdf of C=ln( Z ) is simply the convolution product of the pdfs of A and B.
Once you know the pdf of C it remains simply to compute the one of exp( C ).

See worksheet TRICK.mw for more details.

No solution...

Answered: sand15 1182

September 26 2025

3 3

Unless Maple is wrong or if rewritting the system in a more concise form removes some simplifications or factorizations, it doesn't seem a solution exists:

>

K2 := -(Pn-Pr)/(1-delta)+(-beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+Pr)/delta+(Pr-w-Cr)*beta*(-2+2*delta)*tau*upsilon/(((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))*delta) = 0:

>

K3 := 1-(Pn-Pr)/(1-delta)-(Pn-Cn)/(1-delta)+(Pr-w-Cr)*(1/(1-delta)-(-beta*(Cr*i2-Cr*tau)*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon*Cr/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))^2)/delta)+Ce*rho0/(1-delta) = 0:

>

Error, (in LinearAlgebra:-Determinant) time expired

A good strategy almost always consits in rewritting the equations in a more concise form

>

rewrite := proc(e, v)
  local f1, f2, C, M, f3:
  # assuming that the roots of the numerator do not cancel out the denominator:
  f1 := numer(lhs(e));

  C := [coeffs(f1, [Pn, Pr, w], 'M')];
  M := [M];
  f2 := [seq([seq(degree(t, i), i in [Pn, Pr, w])], t in M)];
  f3 := add(v[f2[i]]*M[i], i=1..nops(M));
  return f3=0, [seq(v[f2[i]] = C[i], i=1..nops(M)) ];
end proc:

newK2, rw2 := rewrite(K2, U):
newK3, rw4 := rewrite(K3, V):
newK4, rw4 := rewrite(K4, W):

Example: K2 rewritten and the corresponding values of the U coefficients:

>	newK2; print(): print~(rw2):

U[[0, 0, 0]] = -Cr*beta*delta^2*i2*upsilon-Cr*beta*delta^2*tau*upsilon-2*beta*delta^2*i2*tau*upsilon+2*beta*delta^2*tau^2*upsilon+2*Cr*beta*delta*i2*upsilon+2*Cr*beta*delta*tau*upsilon+4*beta*delta*i2*tau*upsilon-4*beta*delta*tau^2*upsilon-Cr*beta*i2*upsilon-Cr*beta*tau*upsilon-2*beta*i2*tau*upsilon+2*beta*tau^2*upsilon

(1)

Here is the rewritten system

>	sys := {newK2, newK3, newK4}: for e in sys do e; print(); end do;

Pn^3*V[[3, 0, 0]]+Pn^2*Pr*V[[2, 1, 0]]+Pn^2*w*V[[2, 0, 1]]+Pn*Pr^2*V[[1, 2, 0]]+Pn*Pr*w*V[[1, 1, 1]]+Pr^3*V[[0, 3, 0]]+Pr^2*w*V[[0, 2, 1]]+Pn^2*V[[2, 0, 0]]+Pn*Pr*V[[1, 1, 0]]+Pn*w*V[[1, 0, 1]]+Pr^2*V[[0, 2, 0]]+Pr*w*V[[0, 1, 1]]+w^2*V[[0, 0, 2]]+Pn*V[[1, 0, 0]]+Pr*V[[0, 1, 0]]+w*V[[0, 0, 1]]+V[[0, 0, 0]] = 0

Pn^3*W[[3, 0, 0]]+Pn^2*Pr*W[[2, 1, 0]]+Pn^2*w*W[[2, 0, 1]]+Pn*Pr^2*W[[1, 2, 0]]+Pn*Pr*w*W[[1, 1, 1]]+Pr^3*W[[0, 3, 0]]+Pr^2*w*W[[0, 2, 1]]+Pn^2*W[[2, 0, 0]]+Pn*Pr*W[[1, 1, 0]]+Pn*w*W[[1, 0, 1]]+Pr^2*W[[0, 2, 0]]+Pr*w*W[[0, 1, 1]]+w^2*W[[0, 0, 2]]+Pn*W[[1, 0, 0]]+Pr*W[[0, 1, 0]]+w*W[[0, 0, 1]]+W[[0, 0, 0]] = 0

(2)

newK2 being linear wrt w just solve it wrt w:

>	infolevel[solve] := 4: wsol_1 := solve(newK2, w);

Main: Entering solver with 1 equation in 1 variable

Transformer: solving the uncoupled linear subsystem in w

Main: solving successful - now forming solutions
Main: Exiting solver returning 1 solution

-(Pn^2*U[[2, 0, 0]]+Pn*Pr*U[[1, 1, 0]]+Pr^2*U[[0, 2, 0]]+Pn*U[[1, 0, 0]]+Pr*U[[0, 1, 0]]+U[[0, 0, 0]])/U[[0, 0, 1]]

(3)

Form a 2-by-2 subsystem while plugging the solution above into {newK3 and newK4}

>	sys34 := eval({newK3, newK4}, w=wsol_1): newnewK3, newrw3 := rewrite(sys34[1], X): newnewK4, newrw4 := rewrite(sys34[2], Y):

>	newsys := {newnewK3, newnewK4}: for e in newsys do e; print(); end do;

Pn^4*X[[4, 0, 0]]+Pn^3*Pr*X[[3, 1, 0]]+Pn^2*Pr^2*X[[2, 2, 0]]+Pn*Pr^3*X[[1, 3, 0]]+Pr^4*X[[0, 4, 0]]+Pn^3*X[[3, 0, 0]]+Pn^2*Pr*X[[2, 1, 0]]+Pn*Pr^2*X[[1, 2, 0]]+Pr^3*X[[0, 3, 0]]+Pn^2*X[[2, 0, 0]]+Pn*Pr*X[[1, 1, 0]]+Pr^2*X[[0, 2, 0]]+Pn*X[[1, 0, 0]]+Pr*X[[0, 1, 0]]+X[[0, 0, 0]] = 0

Pn^4*Y[[4, 0, 0]]+Pn^3*Pr*Y[[3, 1, 0]]+Pn^2*Pr^2*Y[[2, 2, 0]]+Pn*Pr^3*Y[[1, 3, 0]]+Pr^4*Y[[0, 4, 0]]+Pn^3*Y[[3, 0, 0]]+Pn^2*Pr*Y[[2, 1, 0]]+Pn*Pr^2*Y[[1, 2, 0]]+Pr^3*Y[[0, 3, 0]]+Pn^2*Y[[2, 0, 0]]+Pn*Pr*Y[[1, 1, 0]]+Pr^2*Y[[0, 2, 0]]+Pn*Y[[1, 0, 0]]+Pr*Y[[0, 1, 0]]+Y[[0, 0, 0]] = 0

(4)

Note that the lhs of each equations are complete polynomials of total degree 4 in the two indeterminates Pn and Pr

>	infolevel[solve] := 4: sol := timelimit(20, solve(newsys, {Pn, Pr}));

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 system split into 1 parts under preprocessing
Main: trying resultant methods
PseudoResultant: 4093009 [2 7000140160 Pn] 2 2 649 0 3 0

PseudoResultant: factoring all equations length = 649
PseudoResultant: factoring done length = 649
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512

LinearEq: subresultant determinant of length 92449 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512
LinearEq: subresultant determinant of length 92449 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 322 and 322

Error, (in convert/sqrfree/sqrfree) time expired

>	timelimit(20, SolveTools:-PolynomialSystem(lhs~(newsys), {Pn, Pr}));

Main: polynomial system split into 1 parts under preprocessing

Main: trying resultant methods
PseudoResultant: 4445613 [2 7000142560 Pn] 2 2 809 0 3 0
PseudoResultant: factoring all equations length = 809

PseudoResultant: factoring done length = 809
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850
LinearEq: subresultant determinant of length 127039 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850

LinearEq: subresultant determinant of length 127039 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 402 and 402
SolutionsLost: setting solutions lost flag
PseudoResultant: 0 solutions found, now doing backsubstitution

SolutionsLost: setting solutions lost flag

Seems to be a dead line

Download Q_solve_sand15.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

1182 Reputation

15 Badges

MaplePrimes Activity

These are answers submitted by sand15

The smallest value of Em is......

Huge errors...

One way...

One way...

Another way...

One way...

value...

Something like that?...

Something like this?...

Explanations...

KKT doesn't seem necessary here...

elementwise comparisons...

Suggestion...

Simply...

No solution...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again