vv

13837 Reputation

20 Badges

9 years, 318 days

MaplePrimes Activity


These are answers submitted by vv

It's possible in any dimension using quadratic forms. Here is an example when the quadric has a center.

[Your question is again more about maths than programming].

restart;

with(LinearAlgebra):

n:=3;
f := 7*x[1]^2 + 4*x[1]*x[2] - 4*x[1]*x[3] + 4*x[2]^2 - 2*x[2]*x[3] + 4*x[3]^2 - 4*x[1] + 5*x[2] + 4*x[3] - 18

3

 

7*x[1]^2+4*x[1]*x[2]-4*x[1]*x[3]+4*x[2]^2-2*x[2]*x[3]+4*x[3]^2-4*x[1]+5*x[2]+4*x[3]-18

(1)

plots:-implicitplot3d(f, x[1]=-5..5, x[2]=-5..5, x[3]=-5..5, style=surface, numpoints=64000, scaling=constrained);

 

A:=VectorCalculus:-Hessian(1/2*f,[seq(x[i],i=1..n)]);

Matrix(3, 3, {(1, 1) = 7, (1, 2) = 2, (1, 3) = -2, (2, 1) = 2, (2, 2) = 4, (2, 3) = -1, (3, 1) = -2, (3, 2) = -1, (3, 3) = 4})

(2)

b:=eval(Vector( [ seq(diff(f,x[i]),i=1..n)]), [seq(x[i]=0,i=1..n)]);

Vector(3, {(1) = -4, (2) = 5, (3) = 4})

(3)

c:=eval(f,[seq(x[i]=0,i=1..n)]);

-18

(4)

X:=Vector([seq(x[i],i=1..n)]);

Vector(3, {(1) = x[1], (2) = x[2], (3) = x[3]})

(5)

solve([ seq(diff(f,x[i]),i=1..n)],{seq(x[i],i=1..n)}); # the center

{x[1] = 11/27, x[2] = -26/27, x[3] = -29/54}

X0:= Vector[column]( eval([seq(x[i],i=1..n)],%) );

Vector(3, {(1) = 11/27, (2) = -26/27, (3) = -29/54})

(6)

J,Q := JordanForm(A,output=['J','Q']);

J, Q := Matrix(3, 3, {(1, 1) = 3, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 9, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 3}), Matrix(3, 3, {(1, 1) = 5/6, (1, 2) = 2/3, (1, 3) = 1/2, (2, 1) = -1/3, (2, 2) = 1/3, (2, 3) = 0, (3, 1) = 4/3, (3, 2) = -1/3, (3, 3) = 1})

T:=Matrix(GramSchmidt([seq( Q[..,j],j=1..n)],normalized)); # T is orthogonal

Matrix(3, 3, {(1, 1) = (5/93)*sqrt(93), (1, 2) = (1/3)*sqrt(6), (1, 3) = -(1/31)*sqrt(62), (2, 1) = -(2/93)*sqrt(93), (2, 2) = (1/6)*sqrt(6), (2, 3) = (7/62)*sqrt(62), (3, 1) = (8/93)*sqrt(93), (3, 2) = -(1/6)*sqrt(6), (3, 3) = (3/62)*sqrt(62)})

(7)

fnew:=simplify( (T.X+X0)^+. A. (T.X+X0) + b.(T.X+X0) + c ); # ==> ellipsoid

-602/27+3*x[1]^2+9*x[2]^2+3*x[3]^2

(8)


 

Download quadric.mw

So, you must help Maple a bit.

restart;

f:=eval(1/(x+2)^4*(-(x+2)^5*RootOf(-2+(x^5*C[1]^5+10*x^4*C[1]^5+40*x^3*C[1]^5+80*x^2*C[1]^5+80*x*C[1]^5+32*C[1]^5)*_Z^25+(5*x^5*C[1]^5+50*x^4*C[1]^5+200*x^3*C[1]^5+400*x^2*C[1]^5+400*x*C[1]^5+160*C[1]^5)*_Z^20)^20*C[1]^5+1)/RootOf(-2+(x^5*C[1]^5+10*x^4*C[1]^5+40*x^3*C[1]^5+80*x^2*C[1]^5+80*x*C[1]^5+32*C[1]^5)*_Z^25+(5*x^5*C[1]^5+50*x^4*C[1]^5+200*x^3*C[1]^5+400*x^2*C[1]^5+400*x*C[1]^5+160*C[1]^5)*_Z^20)^20/C[1]^5, C[1]=t);

(-(x+2)^5*RootOf(-2+(t^5*x^5+10*t^5*x^4+40*t^5*x^3+80*t^5*x^2+80*t^5*x+32*t^5)*_Z^25+(5*t^5*x^5+50*t^5*x^4+200*t^5*x^3+400*t^5*x^2+400*t^5*x+160*t^5)*_Z^20)^20*t^5+1)/((x+2)^4*RootOf(-2+(t^5*x^5+10*t^5*x^4+40*t^5*x^3+80*t^5*x^2+80*t^5*x+32*t^5)*_Z^25+(5*t^5*x^5+50*t^5*x^4+200*t^5*x^3+400*t^5*x^2+400*t^5*x+160*t^5)*_Z^20)^20*t^5)

(1)

r:=indets(f,RootOf)[];

RootOf(-2+(t^5*x^5+10*t^5*x^4+40*t^5*x^3+80*t^5*x^2+80*t^5*x+32*t^5)*_Z^25+(5*t^5*x^5+50*t^5*x^4+200*t^5*x^3+400*t^5*x^2+400*t^5*x+160*t^5)*_Z^20)

(2)

ra:=asympt(r, t, 2);

RootOf(-2+(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160)*_Z^20)*(1/t)^(1/4)-(1/250)*(1/t)^(3/2)/(RootOf(-2+(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160)*_Z^20)^14*(x+2)^5)+O(1/t^2)

(3)

ra:=convert(ra,radical);

2^(1/20)*(1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(1/20)*(1/t)^(1/4)-(1/500)*2^(3/10)*(1/t)^(3/2)/((1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(7/10)*(x+2)^5)+O(1/t^2)

(4)

ff:=eval(f, r=ra);

(-(x+2)^5*(2^(1/20)*(1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(1/20)*(1/t)^(1/4)-(1/500)*2^(3/10)*(1/t)^(3/2)/((1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(7/10)*(x+2)^5)+O(1/t^2))^20*t^5+1)/((x+2)^4*(2^(1/20)*(1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(1/20)*(1/t)^(1/4)-(1/500)*2^(3/10)*(1/t)^(3/2)/((1/(5*x^5+50*x^4+200*x^3+400*x^2+400*x+160))^(7/10)*(x+2)^5)+O(1/t^2))^20*t^5)

(5)

limit(ff, t=infinity);

(3/2)*x+3

(6)

evalf(eval(f=%, [x=2,t=10000])); # check

6.000002820 = 6.

(7)

 

 


 

ineqz:=[abs(z - 1) > 1, abs(z + 1) > 1, Im(z) > 0];

[1 < abs(z-1), 1 < abs(z+1), 0 < Im(z)]

(1)

eval(ineqz, z=x+I*y):

ineqxy:=simplify(%) assuming real;

[1 < (x^2+y^2-2*x+1)^(1/2), 1 < (x^2+y^2+2*x+1)^(1/2), 0 < y]

(2)

plots:-inequal(ineqxy, x=-5..5, y=-5..5);

 


Download ineqz.mw


 

 

 

 

It is always a good idea to simplify the expressions.

restart;

f:=a^2*(2*a^2*p^3-a^2*((p^2+1)^2*(a^2-1))^(1/2)-((p^2+1)^2*(a^2-1))^(1/2)*p^2+2*p*a^2-2*p^3+((p^2+1)^2*(a^2-1))^(1/2)-2*p)/((p^2+1)^2*(a^2-1))^(1/2)/(p^3-((p^2+1)^2*(a^2-1))^(1/2)+p)/(a^2-p^2-1):
# int(f,p); # division by 0

F:=simplify(eval(f, a^2=1+b^2)) assuming b>0,p::real;

(b^2+1)/((b+p)*(p^2+1))

(1)

intf:=eval(int(F, p), b=sqrt(a^2-1));;

-(1/2)*ln(p^2+1)*(a^2-1)/a^2-(1/2)*ln(p^2+1)/a^2+(a^2-1)^(3/2)*arctan(p)/a^2+(a^2-1)^(1/2)*arctan(p)/a^2+ln((a^2-1)^(1/2)+p)*(a^2-1)/a^2+ln((a^2-1)^(1/2)+p)/a^2

(2)

simplify(diff(intf, p) - f) assuming a>0, p::real;  # check;

0

(3)


 

with(plots):

cylx := y^2 + z^2 - 1;
cylz := x^2 + y^2 - 1;

y^2+z^2-1

 

x^2+y^2-1

(1)

cyl:=implicitplot3d([cylx,cylz], x=-2..2, y=-2..2, z=-2..2,
scaling=constrained, style=surface, grid=[40, 40, 40], labels=[x, y, z]):

s:=solve([cylx,cylz],[x,y,z], explicit)

[[x = z, y = (-z^2+1)^(1/2), z = z], [x = z, y = -(-z^2+1)^(1/2), z = z], [x = -z, y = (-z^2+1)^(1/2), z = z], [x = -z, y = -(-z^2+1)^(1/2), z = z]]

(2)

a:=1.05:
ell:=spacecurve([seq(a*rhs~(s[i]),i=1..4)], z=-1..1, color=[red$2,yellow$2], thickness=4):

display(cyl,ell);

 

 

Download cyl.mw

Here is a procedure for Conic classification (general case). It was extracted (& adapted) from a MathApp.
It will not be useful to you, unless you read about this subject.

P:=proc(a,b,c,d,e,f)  #  See ?MathApps,ConicSections
    local equ,disc,Delta, Type, Point, deg, equ2;
    equ := sort(a*x^2+b*x*y+c*y^2+d*x+e*y+f);
    Delta := 4*a*c*f-a*e^2+b*e*d-b^2*f-c*d^2;
    disc := b^2-4*a*c;
    if Delta <> 0 then

        if disc < 0 then
            if c*Delta > 0 then
                Type := "This equation has no real solutions.";
            elif b = 0 and a = c then
                Type := "A circle.";
            else
                Type := "An ellipse.";
            end if;
        elif disc = 0 then
            Type := "A parabola.";
        else
            Type := "A hyperbola.";
        end if;

    # From here on, Delta = 0 : we have the degenerate cases

    elif disc > 0 then #OK
        Type := "Two intersecting lines.";
    elif disc < 0 then #OK
        Type := "A single point.";
        Point := [(2*c*d-b*e)/disc, (2*a*e-b*d)/disc];
    elif equ = 0 then
        Type := "Every point is a solution.";
    else
        deg := degree(equ, [x,y]);
        if deg = 0 then
            Type := "This equation has no solutions.";
        elif deg = 1 then
            Type := "A line.";
            equ2 := equ;
        else
            # a,b,c not all 0
            if c = 0 then
                # c=b=e=0
                disc := d^2-4*a*f;
            else
                # if a = 0 then                
                    # a=b=d=0
                    disc := e^2-4*c*f;
                # else
                    # both discs above have same sign,
                    # so either one is sufficient
                # end if;
            end if;
            if disc > 0 then
                Type := "Two parallel lines.";
            elif disc = 0 then
                Type := "A line.";
                equ2 := select(has, factor(equ)*t, [x,y]);
                if type(equ2, `^`) and op(2, equ2) = 2 then
                    equ2 := op(1, equ2);
                end if;
            else
                Type := "This equation has no solutions.";
            end if;
        end if;
    end if;
    equ=0,'Type'=Type;
end proc:
P(1,2,3,4,5,-6);
plots:-implicitplot(%[1], x=-10..10,y=-10..10, scaling=constrained);
        x^2 + 2*x*y + 3*y^2 + 4*x + 5*y - 6 = 0, Type = "An ellipse."

Why do you say that the first ODE y' = p1/p2 is exact? It became exact, after you have multiplied it with the integrating factor p2.
Most ODEs are in this situation. Should all be declared exact?

## The simplest way is to solve the equation in x
X:=fsolve(1/x^3-sin(12*x)=0, x=1.2..2, maxsols=10); # there are 3 solutions

     X := 1.266058342, 1.591679996, 1.818677859

Y := seq(1/x^3, x=X);

     Y := 0.4927638533, 0.2479891760, 0.1662389018

 

## As a system
s:={fsolve([y-sin(12*x)=0,x^3*y=1], {x,y}, {x=1.2..2})}; # first root

     s := {{x = 1.818677859, y = 0.1662389017}}

fsolve([y-sin(12*x)=0,x^3*y=1], {x,y}, {x=1.2..2}, avoid=s): # repeat if you want
s:=s union {%};

     s := {{x = 1.266058343, y = 0.4927638522}, {x = 1.818677859, y = 0.1662389017}}

N := 100:
for a to N do
for b from a to N do
for c from b to N do
  s:=[isolve(t^3-(a+b+c)*t^2+(a*b+a*c+b*c)*t-2*a*b*c)];
  if nops(s)>1 then print(rhs~(op~(s)), [a,b,c]) fi;
od od od:
                     [3, 5, 16], [1, 8, 15]
                    [7, 8, 30], [2, 15, 28]
                    [6, 10, 32], [2, 16, 30]
                    [5, 16, 42], [2, 21, 40]
                    [9, 15, 48], [3, 24, 45]
                    [8, 21, 45], [3, 35, 36]
                    [7, 33, 80], [3, 40, 77]
                   [14, 16, 60], [4, 30, 56]
                   [12, 20, 64], [4, 32, 60]
                   [11, 24, 70], [4, 35, 66]
                   [10, 32, 84], [4, 42, 80]
                   [15, 25, 80], [5, 40, 75]
                   [13, 35, 96], [5, 48, 91]
                   [21, 24, 90], [6, 45, 84]
                   [18, 30, 96], [6, 48, 90]
                   [16, 42, 90], [6, 70, 72]
                   [18, 55, 112], [7, 88, 90]
 

You have a restart after with(LinearAlgebra).

So, P[0], P[1] are given and you need P[i] for i>1 in terms of L[i,j].
We may suppose P[0] = [0,0] (the origin) and P[1] = [a,0], a>0 (actually a=L[0,1]).
Only L[k,0], L[k,1] will be needed (k>1).

There are 2 solutions (symmetric wrt the x axis), so we may take the y-coordonate of P[k] be positive.

solve( [ x^2+y^2 = L[k,0]^2, (x-a)^2+y^2 = L[k,1]^2 ], {x,y}, explicit ):
select( u -> (sign(eval(y,u))=1), %):
P[k]=eval([x,y],%);

 

 

Use
G1:=RelabelVertices(G,[1,2,3,4,5,6]);

In Maple, a polynomial cannot have symbolic exponents (powers). For example, x^m - x^2 - 2 is not a polynomial, unless m is a (constant) nonnegative integer. This is normal; try for example to compute by hand the quotient and the remainder for the polynomials  x^m - x^2 - 2  and x^3 + 2*x + 3; such operations are essential in the Groebner package.

In conclusion your problem cannot be solved in a CAS (but you may solve it providing integer value(s) for m).

Maple needs a little help.

restart;
eq := x^2+floor(x)-10:
seq( solve({eval(eq, floor(x)=n), n <= x, x < n+1}), n=-5..5);

        {x = -sqrt(14)}, {x = 2*sqrt(2)}

First 38 39 40 41 42 43 44 Last Page 40 of 120