Pinetree

15 Reputation

3 Badges

2 years, 337 days

MaplePrimes Activity


These are questions asked by Pinetree

Quadratic equation ax^2 + bx + c = 0. Coeff's long and ugly expressions in their own right, but I know that a is nonzero, that there is at least one real solution, and that I want the largest.

So, question: How do I get Maple to return the -b/2a PLUS sqrt((b/2a)^2-c/a) solution?

I have several functional equations in equally many unknown functions of at least two variables, plus parameters.  ("collect" works just for single equations, right?)

I know that for certain parameter ranges, all the functions involved will be quadratic, and I know some coefficients are zero.  That gives me some  coefficients to determine.  I want to

  1. specify the functional equations [done in a very primitive low-tech way in the attachment, using atomic variables rather than indices ... have I done correctly?!?] 
  2. get Maple to collect coefficients (the K's and the L's in the attachment; the variables are (y,z))
  3. get Maple to state an equation system these coefficients have to satisfy (these will unfortunately be coupled quadratics)
  4. get Maple to solve that equation system if possible, and if not: to tell me when (= for what parameter values, parameters being the "remaining letters" in the attachment) I have specified enough coefficients
  5. in case of a solution, get Maple to tell me which coefficients are real and positive (for those that are solution of quadratic eq's: whether a positive solution exists)

Phew. I am still a complete newbie. Edit: Attachment link: STcoeff2match.mw where the equations themselves are EQ0, EQ1 and EQ2 at the bottom. Copying and pasting them, they look like this (download STcoeff2pastedEQs.mw)

0 = -r__0*(K__011*y^2+K__022*z^2-K__012*(y-L__1)*(z-L__2)-K__01*(y-L__1)+K__02*(z-L__2))+(-2*K__011*y+m__1+K__012*(z-L__2)+K__01)*((2/3)*c__1*y-(4/3)*K__11*y+(2/3)*`K__12 `*(z-L__2)+(20/9)*K__011*y-(10/9)*K__012*(z-L__2)-(10/9)*K__01-(10/9)*m__1-(1/3)*c__2*z+(2/3)*K__22*z-(1/3)*`K__21 `*(y-L__1)-(16/9)*K__022*z+(8/9)*K__012*(y-L__1)-(8/9)*K__02+(8/9)*m__2)+(-2*K__022*z+m__2+K__012*(y-L__1)-K__02)*((2/3)*c__2*z-(4/3)*K__22*z+(2/3)*`K__21 `*(y-L__1)+(20/9)*K__022*z-(10/9)*K__012*(y-L__1)+(10/9)*K__02-(10/9)*m__2-(1/3)*c__1*y+(2/3)*K__11*y-(1/3)*`K__12 `*(z-L__2)-(16/9)*K__011*y+(8/9)*K__012*(z-L__2)+(8/9)*K__01+(8/9)*m__1)+(-(4/3)*K__011*y+(2/3)*K__022*z+(2/3)*K__012*(z-L__2)-(1/3)*K__012*(y-L__1)-(1/3)*m__2+(2/3)*m__1+(1/3)*K__02+(2/3)*K__01)^2+((2/3)*K__011*y-(4/3)*K__022*z-(1/3)*K__012*(z-L__2)+(2/3)*K__012*(y-L__1)+(2/3)*m__2-(1/3)*m__1-(2/3)*K__02-(1/3)*K__01)^2:

``

0 = -r__1*(K__11*y^2-`K__12 `*y*(z-L__2))+`K__12 `*y*((2/3)*c__2*z-(4/3)*K__22*z+(2/3)*`K__21 `*(y-L__1)+(20/9)*K__022*z-(10/9)*K__012*(y-L__1)+(10/9)*K__02-(10/9)*m__2-(1/3)*c__1*y+(2/3)*K__11*y-(1/3)*`K__12 `*(z-L__2)-(16/9)*K__011*y+(8/9)*K__012*(z-L__2)+(8/9)*K__01+(8/9)*m__1)+((2/3)*c__1*y-(4/3)*K__11*y+(2/3)*`K__12 `*(z-L__2)-(1/3)*c__2*z+(2/3)*K__22*z-(1/3)*`K__21 `*(y-L__1)-(10/9)*K__022*z+(5/9)*K__012*(y-L__1)-(5/9)*K__02+(5/9)*m__2+(8/9)*K__011*y-(4/9)*K__012*(z-L__2)-(4/9)*K__01-(4/9)*m__1)^2:

``

0 = -r__2*(K__22*z^2-`K__21 `*(y-L__1)*z)+`K__21 `*z*((2/3)*c__1*y-(4/3)*K__11*y+(2/3)*`K__12 `*(z-L__2)+(20/9)*K__011*y-(10/9)*K__012*(z-L__2)-(10/9)*K__01-(10/9)*m__1-(1/3)*c__2*z+(2/3)*K__22*z-(1/3)*`K__21 `*(y-L__1)-(16/9)*K__022*z+(8/9)*K__012*(y-L__1)-(8/9)*K__02+(8/9)*m__2)+((2/3)*c__2*z-(4/3)*K__22*z+(2/3)*`K__21 `*(y-L__1)-(1/3)*c__1*y+(2/3)*K__11*y-(1/3)*`K__12 `*(z-L__2)-(10/9)*K__011*y+(5/9)*K__012*(z-L__2)+(5/9)*K__01+(5/9)*m__1+(8/9)*K__022*z-(4/9)*K__012*(y-L__1)+(4/9)*K__02-(4/9)*m__2)^2:

``

 

 

Apologies for possible double post, it seemingly locked up upon trying to post the first time. 

So before I get as far as to ask for how to get a certain PDE system solved and plotted, I tried to fiddle a little bit around with linear algebra. 

  • First, the file SystemGoesWrong.mw . I have issues with declaring (same result if I remove the with(VectorCalculus)); as far as I can see, EQ0 and EQ00 should be the same, except that I have summed the vector in one of them. And the first that "works", is wrong: it returns a scaling of a vector. How come?
  • But then I copy everything from the heading and down into a worksheet where I was already fighting some linear algebra things (can someone please explain?): SystemDeclaresBut.mw 
    Then EQ0 and EQ00 declare just fine! What is the issue?
  • How do I get Maple to list the equations in "compact" form with vector-valued functions so that I can read and debug?  The actual PDE system I want to solve (numerically, of course), looks as follows: DE4Maple.pdf 
    That was also the reason why I tried to declare procedures (coordinate-wise maximum ...), but I guess that questions on how to extract a solution and plot it in a particular way will be its own posting after I have learned how to declare it.

 

The contents of the first file:


 

restart

# Since I do not have any idea of how to get vectors nicely, ...

... I replace U0, U1, U2 by u,v,w and use difftables U,V,W.  And y1, y2 replaced by y,z.

with(PDEtools):

declare(u(y, z), v(y, z), w(y, z), q1(y, z), r1(y, z), s1(y, z), q2(y, z), r2(y, z), s2(y, z), F1(y, z), F2(y, z)):

u(y, z)*`will now be displayed as`*u

 

v(y, z)*`will now be displayed as`*v

 

w(y, z)*`will now be displayed as`*w

 

q1(y, z)*`will now be displayed as`*q1

 

r1(y, z)*`will now be displayed as`*r1

 

s1(y, z)*`will now be displayed as`*s1

 

q2(y, z)*`will now be displayed as`*q2

 

r2(y, z)*`will now be displayed as`*r2

 

s2(y, z)*`will now be displayed as`*s2

 

F1(y, z)*`will now be displayed as`*F1

 

F2(y, z)*`will now be displayed as`*F2

(1)

``

M := `<|>`(`<,>`(2, -1), `<,>`(-1, 2))

M := Matrix(2, 2, {(1, 1) = 2, (1, 2) = -1, (2, 1) = -1, (2, 2) = 2})

(2)

1/M

Matrix([[2/3, 1/3], [1/3, 2/3]])

(3)

 

# The system

 

EQ0 := VectorCalculus:-`+`(VectorCalculus:-`+`(Typesetting:-delayDotProduct(`<|>`(F1, F2), VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))), VectorCalculus:-`-`(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/4, VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))^%T), VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))))), VectorCalculus:-`-`(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/2, `<|>`(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(2, F1), U[y]), VectorCalculus:-`-`(s1)), VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(2, F2), V[z]), VectorCalculus:-`-`(s2)))), 1/M), `<,>`(q1, q2))))

EQ0 := -((1/6)*q1+(1/12)*q2-(1/4)*(diff(u(y, z), y)))*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))-((1/12)*q1+(1/6)*q2-(1/4)*(diff(v(y, z), z)))*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))+(Vector(1, {(1) = F1*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))+F2*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))}, attributes = [coords = cartesian]))+(Vector(1, {(1) = -((2/3)*F1+(1/3)*(diff(u(y, z), y))-(1/3)*s1+(1/3)*F2+(1/6)*(diff(v(y, z), z))-(1/6)*s2)*q1-((1/3)*F1+(1/6)*(diff(u(y, z), y))-(1/6)*s1+(2/3)*F2+(1/3)*(diff(v(y, z), z))-(1/3)*s2)*q2}, attributes = [coords = cartesian]))

(4)

EQ00 := `<|>`(F1, F2).(1/M.`<,>`(q1, q2)-`<,>`(U[y], V[z]))-(1/4)*LinearAlgebra:-Transpose(1/M.`<,>`(q1, q2)-`<,>`(U[y], V[z])).(1/M.`<,>`(q1, q2)-`<,>`(U[y], V[z]))-1/2*(2*`<|>`(F1, F2)+`<|>`(U[y], V[z])-`<|>`(s1, s2)).(1/M).`<,>`(q1, q2)

Error, (in rtable/Sum) invalid input: dimensions do not match: Matrix(1 .. 1, 1 .. 2) cannot be added to Vector[row](1 .. 2)

 

``

``


 

Download SystemGoesWrong.mw

 

 

 

... and of the second:

 

 


 

``

Over to some linear algebra.  

 

 

restart; with(LinearAlgebra); with(VectorCalculus)

NULL

LinearAlgebra:-Transpose(`<|>`(3, 4)).`<,>`(2, 3)

18

(1)

VectorCalculus:-DotProduct(VectorCalculus:-`<,>`(3, 4), VectorCalculus:-`<,>`(2, 3))

18

(2)

DotProduct(`<,>`(2, 3), `<|>`(3, 4))

Matrix([[6, 8], [9, 12]])

(3)

Trace(Matrix(%id = 18446746888362217830));

18

(4)

                                                                 

define(normsqbyDot, normsqbyDot(y::Vector) = VectorCalculus:-DotProduct(y, y))

showstat(normsqbyDot)


normsqbyDot := proc()
local theArgs, arg, look, me, cf, term;
   1   me := eval(procname,1);
   2   theArgs := args;
   3   look := tablelook(('procname')(theArgs),'[`/POS`(1,normsqbyDot,1), `/BIND`(1,1,`/y1`::VectorCalculus:-Vector), `/PATTERN`(`/y1`^2)]');
   4   if look <> FAIL then
   5     eval(look,`/FUNCNAME` = procname)
       else
   6     ('procname')(theArgs)
       end if
end proc

 

define, "%1 is assigned", normsqbyMatrixProduct

showstat(normsqbyMatrixProduct)


normsqbyMatrixProduct := proc()
local theArgs, arg, look, me, cf, term;
   1   me := eval(procname,1);
   2   theArgs := args;
   3   look := tablelook(('procname')(theArgs),'[`/POS`(1,normsqbyMatrixProduct,1), `/BIND`(1,1,`/y1`::Matrix), `/PATTERN`(`/y1`^2)]');
   4   if look <> FAIL then
   5     eval(look,`/FUNCNAME` = procname)
       else
   6     ('procname')(theArgs)
       end if
end proc

 

normsqbyDot(`<,>`(2, 3))

Error, (in rtable/Power) exponentiation operation not defined for Vectors

 

normsqbyMatrixProduct(VectorCalculus:-`<,>`(2, 3))

normsqbyMatrixProduct(Vector(2, {(1) = 2, (2) = 3}, attributes = [coords = cartesian]))

(5)

convert(n*ormsqbyMatrixProduct(`<,>`(2, 3)), float)

n*ormsqbyMatrixProduct(Vector(2, {(1) = 2, (2) = 3}, attributes = [coords = cartesian]))

(6)

NULL

# Since I do not have any idea of how to get vectors nicely, ...

... I replace U0, U1, U2 by u,v,w and use difftables U,V,W.  And y1, y2 replaced by y,z.

with(PDEtools); -1; with(plots)

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

(7)

declare(u(y, z), v(y, z), w(y, z), q1(y, z), r1(y, z), s1(y, z), q2(y, z), r2(y, z), s2(y, z), F1(y, z), F2(y, z)):

u(y, z)*`will now be displayed as`*u

 

v(y, z)*`will now be displayed as`*v

 

w(y, z)*`will now be displayed as`*w

 

q1(y, z)*`will now be displayed as`*q1

 

r1(y, z)*`will now be displayed as`*r1

 

s1(y, z)*`will now be displayed as`*s1

 

q2(y, z)*`will now be displayed as`*q2

 

r2(y, z)*`will now be displayed as`*r2

 

s2(y, z)*`will now be displayed as`*s2

 

F1(y, z)*`will now be displayed as`*F1

 

F2(y, z)*`will now be displayed as`*F2

(8)

NULL

M := `<|>`(`<,>`(2, -1), `<,>`(-1, 2))

M := Matrix(2, 2, {(1, 1) = 2, (1, 2) = -1, (2, 1) = -1, (2, 2) = 2})

(9)

1/M

Matrix([[2/3, 1/3], [1/3, 2/3]])

(10)

 

# The system

 

EQ0 := Typesetting:-delayDotProduct(`<|>`(F1, F2), Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2))+`-`(`<,>`(U[y], V[z])))+`-`(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/4, (Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2))+`-`(`<,>`(U[y], V[z])))^%T), Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2))+`-`(`<,>`(U[y], V[z]))))+`-`(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/2, `<|>`(VectorCalculus:-`*`(2, F1)+U[y]+`-`(s1), VectorCalculus:-`*`(2, F2)+V[z]+`-`(s2))), 1/M), `<,>`(q1, q2)))

F1*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))+F2*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))-((1/6)*q1+(1/12)*q2-(1/4)*(diff(u(y, z), y)))*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))-((1/12)*q1+(1/6)*q2-(1/4)*(diff(v(y, z), z)))*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))-((2/3)*F1+(1/3)*(diff(u(y, z), y))-(1/3)*s1+(1/3)*F2+(1/6)*(diff(v(y, z), z))-(1/6)*s2)*q1-((1/3)*F1+(1/6)*(diff(u(y, z), y))-(1/6)*s1+(2/3)*F2+(1/3)*(diff(v(y, z), z))-(1/3)*s2)*q2

(11)

EQ00 := VectorCalculus:-`+`(VectorCalculus:-`+`(Typesetting:-delayDotProduct(`<|>`(F1, F2), VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))), VectorCalculus:-`-`(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/4, VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))^%T), VectorCalculus:-`+`(Typesetting:-delayDotProduct(1/M, `<,>`(q1, q2)), VectorCalculus:-`-`(`<,>`(U[y], V[z])))))), VectorCalculus:-`-`(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(VectorCalculus:-`*`(1/2, VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(2, `<|>`(F1, F2)), `<|>`(U[y], V[z])), VectorCalculus:-`-`(`<|>`(s1, s2)))), 1/M), `<,>`(q1, q2))))

F1*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))+F2*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))-((1/6)*q1+(1/12)*q2-(1/4)*(diff(u(y, z), y)))*((2/3)*q1+(1/3)*q2-(diff(u(y, z), y)))-((1/12)*q1+(1/6)*q2-(1/4)*(diff(v(y, z), z)))*((1/3)*q1+(2/3)*q2-(diff(v(y, z), z)))-((2/3)*F1+(1/3)*(diff(u(y, z), y))-(1/3)*s1+(1/3)*F2+(1/6)*(diff(v(y, z), z))-(1/6)*s2)*q1-((1/3)*F1+(1/6)*(diff(u(y, z), y))-(1/6)*s1+(2/3)*F2+(1/3)*(diff(v(y, z), z))-(1/3)*s2)*q2

(12)

``

 


 

Download SystemDeclaresBut.mw

 

 

So I encountered a Chini differential equation which is solvable in terms of an implicit function

restart; assume(nu > 2, beta > 0, lambda > 0, delta > 0, y >= 0);
simplify(dsolve(beta+y+(D(h))(y) = nu*((lambda+beta+y)^2-delta*h(y))^(1/2)));

The answer is really too ugly to be pasted here, but I do get a
RootOf [ linear + log(linear) + delta* some integral + 2 * the very same integral]

It took me a quarter of an hour to deciphre that the integrals are the same and that I could get it all down to (delta+2)*integral. Because Maple did not simplify it. 

But if I copy what is inside the RootOf(), and paste it into a simplify() - then it does! So how to force it into doing that in the first place? 

 

Still novice, yes ... 

So I needed a CAS, and I spent a couple of months trying to get a basic understanding of SymPy and various applications starting with "M". We have Maple version 2016.

My "prototype problem" can be solved by hand, and is a system of quadratic Bellman equations, for i=1,2.  I'll return to it below, as I am obviously too clumsy to get even the second-to-simplest max/min working.  Oh, and I can't even insert Maple Math here in the forum, it does not like maximize or minimize.

Let's start easy. I enter
maximize((b-x)x,x) 
which works as I expect. Then already at
maximize((abs(b)-x)x,x>=0) 
I am stuck. Please, sweet Maple, you know that the answer is the same as before, don't you? (In the meantime I have tried to feed it maximize((b-x)x,x=0..1)... )

I was hoping my "proper" problem should be doable. What I really need is a sequence of quadratic optimization problems, where I have a vector x maximizing b'x-x'Ax subject to linear constraints, so it should not be too hard.  The "prototype problem" I needed for starters, is a system where for i=1,2 I have 

v_i(0,y)=0 and inductively v_i(t+1,y)=max{x_i*(b_i-x_1-x_2)+ r_i v_i(t,y-x)}

where b_i and r_i are constants, x=(x_1,x_2), y=(y_1,y_2), and everything is nonnegative - including, the choice variables x_i must be between 0 and min{y_i, b_i-x_1-x_2}. And I want to plot both functions and the x with time as a slider, but ... I don't think I'll ever get that far? 

Sorry for whining. (I know why I dropped out of computer science.)

Page 1 of 1