clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Numeri o lettere vengono 
%% definiti come simboli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = sqrt(2)
y = sqrt(sym(2))
double(y)
z = sqrt(sym(3))
y + z

a=10/3+9/7
b=sym(10)/sym(3)+sym(9)/sym(7)
double(b)

syms a b
c = (a+b)^2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Semplificazione di espressioni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms a b
g = (2*b^2+a*b-a^2)/(a+b)
simplify(g)
simple(g)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fattorizzazione
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms a b
g = (2*b^2+a*b-a^2)/(a+b)
factor(g)

syms x y
factor( x^2-y^2 )
factor( x^3+y^3 )
factor( 2*x^3*y+x^2*y^2+2*x*y^3+y^4 )

factor(sym(272949548788224000))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Raccoglimento
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x y
z = 3*x*y-2*y^2+5*x^2*y-7*x*y^2+11*x^2*y^2
collect(z,x)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sostituire una variabile simbolica
%% con un numero. Il comando subs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms a b c 
g = (a+b)^2-(b+c)^2-(a+c)^2
subs(g, a, 3)

syms n k
bincoef = sym('n!/(k!*(n-k)!)')
subs(bincoef,{n,k},{6,3})

coef = zeros(1,7);
for j=0:6
 coef(j+1)=subs(bincoef,{n,k},{6,j});
end
coef


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Risolvere equazioni e sistemi di equazioni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x
f = x^2 - 7*x + 12
solve(f)

f = 2*x^2 + x + 2
sol = solve(f)
%%%%%
%% Valuta il risultato in termini 
%% numerici, non pił simbolici
%%%%%
vpa(sol)

f1 = x^2+y^2-50
f2 = x*y-7
sol = solve(f1,f2)
isstruct(sol)
sol.x
sol.y

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Soluzione per via simbolica
%% del sistema:
%% x^2+5*y-7*x+7=0
%% y^2+x*y-x^2*y+5=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x y
f1 = x^2+5*y-7*x+7
f2 = y^2+x*y-x^2*y+5
sol = solve(f1,f2)
sol.x
sol.y
vpa(sol.x)
vpa(sol.y)
%%%%
%% Verifica
%% Si sostiusce la soluzione 
%% nel sistema di equazioni
%%%%
for j=1:numel(sol.x)
  simplify([subs(f1,{x,y},{sol.x(j), sol.y(j)})  subs(f2,{x,y},{sol.x(j), sol.y(j)})])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Soluzione per via simbolica
%% del sistema:
%% x^2+5*y-7*x+7=0
%% y^2+x*y-x^2*y+5=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f1 = @(x) [x(1)^2+5*x(2)-7*x(1)+7 x(2)^2+x(1)*x(2)-x(1)^2*x(2)+5]
% A seconda della condizione iniziale
% scelta otteniamo una delle soluzioni
% (ovviamente solo quelle reali)
fsolve(f1,[2.1 1.2])
vpa([sol.x(1) sol.y(1)])
%
fsolve(f1,[5 .1])
vpa([sol.x(2) sol.y(2)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Nello spazio delle condizioni
%% iniziali esistono quindi dei
%% bacini di attrazione
%% Cerchiamo di rappresentarli
%% graficamente (richiede tempo)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warning off;
opts = optimset('Display','off');
n = 21;
x1 = linspace(5,30,n);
x2 = linspace(-20,20,n);
[X1,X2]=meshgrid(x1,x2);
sol = zeros(n);
for r=1:n
    for c=1:n
        out=fsolve(f1,[X1(r,c) X2(r,c)],opts);
        sol(r,c)=out(1);
    end
end
warning on;
figure(20);
imagesc(sol(:,:))

%%%%%%%%%%%%%%%%%%%%%%
%% I coeff. del sistema 
%% sono ora simboli
%% e non valori numerici
%%%%%%%%%%%%%%%%%%%%%%%
syms x y a b
f1 = x^2+y^2-a
f2 = x*y-b
sol = solve(f1,f2,x,y)
sol.x
sol.y
% Verifica
for j=1:numel(sol.x)
  simplify([subs(f1,{x,y},{sol.x(j), sol.y(j)})  subs(f2,{x,y},{sol.x(j), sol.y(j)})])
end

%%%%%%%%%%%%%%%%%%%%%%
%% Sistema impossibile
%%%%%%%%%%%%%%%%%%%%%%%
% Soluzione per via simbolica
syms x y
eq1 = x^2-y^2
eq2 = (x^2-y^2)/(x-y)-sym(1)
sol = solve(eq1,eq2,x,y)
% Attenzione !!! Soluzione errata !
sol.x
sol.y
for j=1:numel(sol.x)
  simplify([subs(eq1,{x,y},{sol.x(j), sol.y(j)})  subs(eq2,{x,y},{sol.x(j), sol.y(j)})])
end

% Soluzione per via numerica
f1 = @(x) [x(1)^2-x(2)^2 (x(1)^2-x(2)^2)/(x(1)-x(2))-1]
f1([sol.x sol.y])
% Attenzione ai warning !
fsolve(f1,[1 1])

%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Curve di livello della
%% funzione di utilitą Ua
%%%%%%%%%%%%%%%%%%%%%%%%%%
syms xa xb k
Ua = log(xa)+5*log(xb)
solve(Ua-k,xa)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Derivazione di funzioni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x y

g = atan(x)
diff(g,x)

f = x*log(x*y)
diff(f,x)
diff(f,y)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Integrazione di funzioni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x y
f = x+1/x+1/x^2+sin(2*x)
int(f)

f = sin(x)/x
int(f)

g = exp(-x^2)
int(g,-Inf,Inf)

f = x*log(x*y)
diff(f,x)
diff(f,y)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Interazione con Maple
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pulisce lo spazio di lavoro di Maple
% e reinizializza il suo kernel
maple restart

% Help di un comando di Maple
mhelp taylor
mhelp taylor/details

% Si richiama il comando taylor di
% Maple con dati parametri
syms x y 
maple('taylor','f(x)','x=0',3)
maple('taylor','f(x)','x=1',3)
maple('taylor',exp(x),'x=0',5)

% Si richiama il comando mtaylor di
% Maple con dati parametri
mhelp mtaylor
maple('mtaylor','f(x,y)','[x=0, y=0]',2)
maple('mtaylor','f(x,y)','[x=0, y=0]',3)
maple('mtaylor',sin(x*y),'[x=0, y=0]',10)
% Attenzione ! Soluzione non corretta.
maple('mtaylor',sin(x*y)/(x*y),'[x=0, y=0]',10)
%%%%%%%%%%%%%%%

% Si richiama il comando sum di
% Maple con dati parametri
syms x k
S = maple('sum',x^k,'k=0..n')
simplify(S)
S = maple('sum',x^k,'k=0..Inf')
S = simplify(maple('sum',x^(-k),'k=0..n'))
S = maple('sum',x^(-k),'k=0..Inf')

% Ci sono funzioni di Maple che non sono
% utilizzabili in Matlab
% Ad esempio plot
% maple('plot',cos(x/2) + sin(2*x),'x = 0..4*Pi')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Soluzione di equazioni differenziali
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dsolve('D2y+y=0')

dsolve('Dy+x^3=0','x')

dsolve('D2y+y-exp(x)=0','x')

dsolve('D3y+y*Dy-x=0','x')

sol = dsolve('x*D2y+(1+x)*Dy=1/(16*x)*y','x')

sol = dsolve('Dy+p(x)*y=q(x)','x')

sol = dsolve('Dx=4*x-3*y','Dy=-2*x-y')
sol.x
sol.y

% Qui si specificano le condizioni al contorno
sol = dsolve('Dx=2*x','Dy=3*x+2*y','x(0)=1','y(0)=1')
sol.x
sol.y

% Tracciamo la soluzione
sol = dsolve('Dx=y','Dy=-0.5*x-y','x(0)=1','y(0)=1')
sol.x
sol.y
tvec = linspace(0,20,300);
X = subs(sol.x,'t',linspace(0,20,200));
Y = subs(sol.y,'t',linspace(0,20,200));
figure(1);
plot(X,Y)

% Tracciamo un fascio di traiettorie
figure(2)
hold on;
axis([0 8 -2.5 2.5])
xx = linspace(0.5,8,100);
for x0=-2:1:2
    cond = ['y(0.5)=',num2str(x0)];
    sol = dsolve('x*Dy+3*y=x*sin(x)',cond,'x');
    yy = subs(sol,'x',xx);
    plot(xx,yy)
    drawnow
end

% Tracciamo un fascio di traiettorie
figure(3)
axis([-1.5 1.5 -1.5 1.5])
hold on;
tt = linspace(0,10,50);
for x0 = -1:.5:1
    for y0 = -1:.5:1
        cond1 = ['x(0)=',num2str(x0)];
        cond2 = ['y(0)=',num2str(y0)];
        sol = dsolve('Dx=y','Dy=-0.5*x-y',cond1,cond2);
        xx = subs(sol.x,'t',tt);
        yy = subs(sol.y,'t',tt);
        plot(xx,yy);
        drawnow
    end
end
hold off

% Integrazione numerica di un'equazione
% differenziale nonlineare con il comando ode45
% Sistema caotico di Lorenz
figure(4);
SIGMA = 10.;
RHO = 28.;
BETA = 8./3.;
ydot = @(t,y) [-BETA*y(1)+y(2)*y(3); -SIGMA*y(2)+SIGMA*y(3); -y(1)*y(2)+RHO*y(2)-y(3)];
y0=[20 5 -5];
[t, traj] = ode45(ydot,[0 50],y0);
plot3(traj(:,1),traj(:,2),traj(:,3),'-')
box on

% I comandi comet e comet3 permettono di 
% ottenere una animazione della traiettoria
figure(5);
comet3(traj(:,1),traj(:,2),traj(:,3))