%Priklad 6.3 i)
syms x y
f=x^4+y^4-x^2-2*x*y-y^2
f =
x^4+y^4-x^2-2*x*y-y^2
fx=diff(f,x)
fx =
4*x^3-2*x-2*y
fy=diff(f,y)
fy =
4*y^3-2*x-2*y
help solve
 SOLVE  Symbolic solution of algebraic equations.
    SOLVE('eqn1','eqn2',...,'eqnN')
    SOLVE('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
    SOLVE('eqn1','eqn2',...,'eqnN','var1','var2',...'varN')

    The eqns are symbolic expressions or strings specifying equations.  The
    vars are symbolic variables or strings specifying the unknown variables.
    SOLVE seeks zeros of the expressions or solutions of the equations.
    If not specified, the unknowns in the system are determined by FINDSYM.
    If no analytical solution is found and the number of equations equals
    the number of dependent variables, a numeric solution is attempted.

    Three different types of output are possible.  For one equation and one
    output, the resulting solution is returned, with multiple solutions to
    a nonlinear equation in a symbolic vector.  For several equations and
    an equal number of outputs, the results are sorted in lexicographic
    order and assigned to the outputs.  For several equations and a single
    output, a structure containing the solutions is returned.

    Examples:

       solve('p*sin(x) = r') chooses 'x' as the unknown and returns

         ans =
         asin(r/p)

       [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0') returns

         x =
         [ 1]
         [ 3]

         y =
         [    1]
         [ -3/2]

       S = solve('x^2*y^2 - 2*x - 1 = 0','x^2 - y^2 - 1 = 0') returns
       the solutions in a structure.

         S =
           x: [8x1 sym]
           y: [8x1 sym]

       [u,v] = solve('a*u^2 + v^2 = 0','u - v = 1') regards 'a' as a
       parameter and solves the two equations for u and v.

       S = solve('a*u^2 + v^2','u - v = 1','a,u') regards 'v' as a
       parameter, solves the two equations, and returns S.a and S.u.

       [a,u,v] = solve('a*u^2 + v^2','u - v = 1','a^2 - 5*a + 6') solves
       the three equations for a, u and v.

    See also <a href="matlab:help dsolve">dsolve</a>.

    Overloaded functions or methods (ones with the same name in other directories)
       <a href="matlab:help sym/solve.m">help sym/solve.m</a>

    Reference page in Help browser
       <a href="matlab:doc solve">doc solve</a>

solve(fx,fy)
ans =
    x: [9x1 sym]
    y: [9x1 sym]
res=solve(fx,fy)
res =
    x: [9x1 sym]
    y: [9x1 sym]
res.x
ans =
                                                0
                                                0
                                                0
                                               -1
                                                1
 -1/4*(1+i*3^(1/2))^(3/2)+1/2*(1+i*3^(1/2))^(1/2)
  1/4*(1+i*3^(1/2))^(3/2)-1/2*(1+i*3^(1/2))^(1/2)
 -1/4*(1-i*3^(1/2))^(3/2)+1/2*(1-i*3^(1/2))^(1/2)
  1/4*(1-i*3^(1/2))^(3/2)-1/2*(1-i*3^(1/2))^(1/2)
res.y
ans =
                        0
                        0
                        0
                       -1
                        1
 -1/2*(1+i*3^(1/2))^(1/2)
  1/2*(1+i*3^(1/2))^(1/2)
 -1/2*(1-i*3^(1/2))^(1/2)
  1/2*(1-i*3^(1/2))^(1/2)
% Ruzna realna reseni jsou body
% P1=[1,1], P2=[-1,-1], P3=[0,0]
% Spocteme Hessovu matici:
fxx=diff(fx,x)
fxx =
12*x^2-2
fyy=diff(fy,y)
fyy =
12*y^2-2
fxy=diff(fx,y)
fxy =
-2
H=[fxx,fxy;fxy,fyy]
H =
[ 12*x^2-2,       -2]
[       -2, 12*y^2-2]
help sym/eig
 EIG    Symbolic eigenvalues and eigenvectors.
    With one output argument, LAMBDA = EIG(A) is a symbolic vector
    containing the eigenvalues of a square symbolic matrix A.

    With two output arguments, [V,D] = EIG(A) returns a matrix V whose
    columns are eigenvectors and a diagonal matrix D containing eigenvalues.
    If the resulting V is the same size as A, then A has a full set of
    linearly independent eigenvectors which satisfy A*V = V*D.

    With three output arguments, [V,D,P] also returns P, a vector of indices
    whose length is the total number of linearly independent eigenvectors,
    so that A*V = V*D(P,P).  If A is n-by-n, then V is n-by-m where n is
    the sum of the algebraic multiplicities and m is the sum of the geometric
    multplicities.

    LAMBDA = EIG(VPA(A)) and [V,D] = EIG(VPA(A)) compute numeric eigenvalues
    and eigenvectors using variable precision arithmetic.  If A does not
    have a full set of eigenvectors, the columns of V will not be linearly
    independent.

    Examples:
       [v,lambda] = eig([a,b,c; b,c,a; c,a,b])

       R = sym(gallery('rosser'));
       eig(R)
       [v,lambda] = eig(R)
       eig(vpa(R))
       [v,lambda] = eig(vpa(R))

       A = sym(gallery(5)) does not have a full set of eigenvectors.
       [v,lambda,p] = eig(A) produces only one eigenvector.

    See also <a href="matlab:help poly">poly</a>, <a href="matlab:help jordan">jordan</a>, <a href="matlab:help svd">svd</a>, <a href="matlab:help vpa">vpa</a>.

lambda=eig(H)
lambda =
 6*x^2-2+6*y^2+2*(9*x^4-18*x^2*y^2+1+9*y^4)^(1/2)
 6*x^2-2+6*y^2-2*(9*x^4-18*x^2*y^2+1+9*y^4)^(1/2)
help subs
 SUBS   Symbolic substitution.
    SUBS(S) replaces all the variables in the symbolic expression S with
    values obtained from the calling function, or the MATLAB workspace.

    SUBS(S,NEW) replaces the free symbolic variable in S with NEW.
    SUBS(S,OLD,NEW) replaces OLD with NEW in the symbolic expression S.
    OLD is a symbolic variable, a string representing a variable name, or
    a string (quoted) expression. NEW is a symbolic or numeric variable
    or expression.

    If OLD and NEW are vectors or arrays of the same size, each element
    of OLD is replaced by the corresponding element of NEW.  If S and OLD
    are scalars and NEW is an array or cell array, the scalars are expanded
    to produce an array result.  If NEW is a cell array of numeric matrices,
    the substitutions are performed elementwise (i.e., subs(x*y,{x,y},{A,B})
    returns A.*B when A and B are numeric).

    If SUBS(S,OLD,NEW) does not change S, then SUBS(S,NEW,OLD) is tried.
    This provides backwards compatibility with previous versions and
    eliminates the need to remember the order of the arguments.
    SUBS(S,OLD,NEW,0) does not switch the arguments if S does not change.

    Examples:
      Single input:
        Suppose a = 980 and C1 = 3 exist in the workspace.
        The statement
           y = dsolve('Dy = -a*y')
        produces
           y = exp(-a*t)*C1
        Then the statement
           subs(y)
        produces
           ans = 3*exp(-980*t)

      Single Substitution:
        subs(a+b,a,4) returns 4+b.

      Multiple Substitutions:
        subs(cos(a)+sin(b),{a,b},[sym('alpha'),2]) or
        subs(cos(a)+sin(b),{a,b},{sym('alpha'),2}) returns
        cos(alpha)+sin(2)

      Scalar Expansion Case:
        subs(exp(a*t),'a',-magic(2)) returns

        [   exp(-t), exp(-3*t)]
        [ exp(-4*t), exp(-2*t)]

      Multiple Scalar Expansion:
        subs(x*y,{x,y},{[0 1;-1 0],[1 -1;-2 1]}) returns
        [  0, -1]
        [  2,  0]

    See also SUBEXPR

    Overloaded functions or methods (ones with the same name in other directories)
       <a href="matlab:help sym/subs.m">help sym/subs.m</a>

    Reference page in Help browser
       <a href="matlab:doc subs">doc subs</a>

subs(lambda,{x,y},[1,1])
ans =
    12
     8
% 2.dif. je pozitivne definitini forma v P1 => P1=lok.min
subs(lambda,{x,y},[-1,-1])
ans =
    12
     8
% Podobne P2
subs(lambda,{x,y},[0,0])
ans =
     0
    -4
% V P3 je forma negativne semidefinitni a o extremu nelze rozhodnout
% Priklad 6.5
syms z
u=x+y^2/(4*x)+z^2/y+2/z
u =
x+1/4*y^2/x+z^2/y+2/z
ux=diff(u,x)
ux =
1-1/4*y^2/x^2
uy=diff(u,y)
uy =
1/2*y/x-z^2/y^2
uz=diff(u,z)
uz =
2*z/y-2/z^2
res=solve(ux,uy,uz)
res =
    x: [8x1 sym]
    y: [8x1 sym]
    z: [8x1 sym]
res.x
ans =
                1/2*i
               -1/2*i
                 -1/2
                  1/2
 (-1/4+1/4*i)*2^(1/2)
  (1/4-1/4*i)*2^(1/2)
 (-1/4-1/4*i)*2^(1/2)
  (1/4+1/4*i)*2^(1/2)
res.y
ans =
                    i
                   -i
                   -1
                    1
  (1/2-1/2*i)*2^(1/2)
 (-1/2+1/2*i)*2^(1/2)
  (1/2+1/2*i)*2^(1/2)
 (-1/2-1/2*i)*2^(1/2)
res.z
ans =
                         -i
                          i
                         -1
                          1
 -1/2*2^(1/2)-1/2*i*2^(1/2)
  1/2*2^(1/2)+1/2*i*2^(1/2)
 -1/2*2^(1/2)+1/2*i*2^(1/2)
  1/2*2^(1/2)-1/2*i*2^(1/2)
% Existuje dve realna reseni, z nichz jen
% bod P=[1/2,1,1] lezi v 1.kvadrantu
uxx=diff(ux,x)
uxx =
1/2*y^2/x^3
uyy=diff(uy,y)
uyy =
1/2/x+2*z^2/y^3
uzz=diff(uz,z)
uzz =
2/y+4/z^3
uxy=diff(ux,y)
uxy =
-1/2*y/x^2
uxz=diff(ux,z)
uxz =
0
uyz=diff(uy,z)
uyz =
-2*z/y^2
H=[uxx,uxy,uxz;uxy,uyy,uyz;uxz,uyz,uzz]
H =
[     1/2*y^2/x^3,      -1/2*y/x^2,               0]
[      -1/2*y/x^2, 1/2/x+2*z^2/y^3,        -2*z/y^2]
[               0,        -2*z/y^2,       2/y+4/z^3]
HvP=subs(H,{x,y,z},[1/2,1,1])
HvP =
     4    -2     0
    -2     3    -2
     0    -2     6
eig(sym(HvP))
ans =
                                                                                                  1/3*(-62+93*sqrt(-1)*3^(1/2))^(1/3)+31/3/(-62+93*sqrt(-1)*3^(1/2))^(1/3)+13/3
 -1/6*(-62+93*sqrt(-1)*3^(1/2))^(1/3)-31/6/(-62+93*sqrt(-1)*3^(1/2))^(1/3)+13/3+1/2*sqrt(-1)*3^(1/2)*(1/3*(-62+93*sqrt(-1)*3^(1/2))^(1/3)-31/3/(-62+93*sqrt(-1)*3^(1/2))^(1/3))
 -1/6*(-62+93*sqrt(-1)*3^(1/2))^(1/3)-31/6/(-62+93*sqrt(-1)*3^(1/2))^(1/3)+13/3-1/2*sqrt(-1)*3^(1/2)*(1/3*(-62+93*sqrt(-1)*3^(1/2))^(1/3)-31/3/(-62+93*sqrt(-1)*3^(1/2))^(1/3))
eval(eig(sym(HvP)))
ans =
   7.2971 - 0.0000i
   0.9161 - 0.0000i
   4.7868 + 0.0000i
% Vsechna vl.c. kladna => poz.def. => lok.min v P

% Priklad 6.6iii)
P=x^2+pi*y^2+(l-4*x-2*pi*y)^2/(12*sqrt(3))
??? Undefined function or variable 'l'.

syms l
P=x^2+pi*y^2+(l-4*x-2*pi*y)^2/(12*sqrt(3))
P =
x^2+pi*y^2+1/36*(l-4*x-2*pi*y)^2*3^(1/2)
Px=diff(P,x)
Px =
2*x-2/9*(l-4*x-2*pi*y)*3^(1/2)
Py=diff(P,y)
Py =
2*pi*y-1/9*(l-4*x-2*pi*y)*3^(1/2)*pi
res=solve(Px,Py)
res =
    x: [1x1 sym]
    y: [1x1 sym]
res.x
ans =
l*(-3*3^(1/2)+pi+4)/(-11+8*pi+pi^2)
res.y
ans =
1/2*l*(-3*3^(1/2)+pi+4)/(-11+8*pi+pi^2)
H1=subs(P,{y}, 0)
H1 =
x^2+1/36*(l-4*x)^2*3^(1/2)
H2=subs(P,x,0)
H2 =
pi*y^2+1/36*(l-2*pi*y)^2*3^(1/2)
H3=subs(P,y,(l-4*x)/(2*pi))
H3 =
x^2+1/4/pi*(l-4*x)^2
H1x=diff(H1,x)
H1x =
2*x-2/9*(l-4*x)*3^(1/2)
H2y=diff(H2,y)
H2y =
2*pi*y-1/9*(l-2*pi*y)*3^(1/2)*pi
H3x=diff(H3,x)
H3x =
2*x-2*(l-4*x)/pi
r1=solve(H1x)
r1 =
3^(1/2)*l/(9+4*3^(1/2))
