%%  CompMath-Vorlesung 13.11.2020
% Demonstration symbolischer Methoden in Matlab, 
% Siehe auch die <http://ewh.ieee.org/r1/ct/sps/PDF/MATLAB/chapter4.pdf gute Einfuehrung> dazu.

%% Lade das symbolische Paket fuer Octave
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if (isOctave)
   setenv PYTHON python
   pkg load symbolic
end

%% Symbolische Funktionen
clear; clc;

disp('Symbolische Funktionen');
%  Deklaration
syms a b c;        % deklariere symbolische Variablen
f = a+b+c          % symbolische Funktion

% Ersetzen von Variablen durch Zahlenwerten oder Ausdruecke
subs(f, a, 3)      % Ersetze a durch 3

subs(f, a, b)      % Ersetze a durch (bereits als symb. def.) Variable b

subs(f, a, 'x')    % Ersetze a durch (noch nicht als symb.) Variable x

subs(f, a, (b+c)^2)  % Ersetze a durch den Term (b+c)^2

subs(f, [a,b,c], [1,2,3])  % Ersetze die Variablen durch die Werte/Terme

%% Manipulation algebraischer Ausdruecke 
%   a^2-2*a*b+b^2
clear; 
%clc;

disp('Manipulation algebraischer Ausdruecke');
syms a b;          % deklariere symb. Variablen
f = a^2-2*a*b+b^2 % symb. Funktion

fe = factor(f)    % faktorisieren

expand(fe)        % ausmultiplizieren

%% simple anwenden
disp('simple anwenden')
% simple(f)         % (bestmoeglich) vereinfachen, zeigt alle Versuche
ff=simplify(f)      % (bestmoeglich) vereinfachen, zeigt alle Versuche
expand(ff)          % ausmultiplizieren


%% Algebraische Gleichungen
% Loese x^2=16
%% Loese nach der einzigen Variablen auf: symbolisch
clear; 
%clc;
disp('Loese nach der einzigen Variablen auf')
% Loese x^2=16:  Variante 1
syms x
solve(x^2 == 16)

% Loese x^2=16:  Variante 2
syms x;
y = x^2-16;
solve(y==0, x)       % Loest  y(x) = 0 

%% Loese symbolisch mit Zusatzannahme
clear;
syms x
assume(x,'clear')                   % ensure that all old assumtions are deleted
f = x^4 - x^3 - 5*x^2 - 7*x - 84
disp('Loese ohne Zusatzannahme')
xsol = solve(f==0)

disp('x soll reell sein')
assume(x,'real')
xsol = solve(f==0)

disp('x soll negativ sein')
assumeAlso(x,'negative')
xsol = solve(f==0)

disp('assumptions(x) :')
assumptions(x)


%% Loese nach der einzigen Variablen auf: numerisch
% Numerisch mit  fzero(), wobei die symbolische Funktion 
% mittels <http://de.mathworks.com/help/symbolic/matlabfunction.html matlabFunction()>
% in eine numerische Funktion konvertiert werden muss.
syms x;
y = x^2-16;

disp('Numerisch mit  fzero():')

x0  = -1;
% 
%f_hand = @(x)subs(y,x);          % Functionshandle aus symbolischer Fkt. erzeugen (alt)
f_hand = matlabFunction(y);       % numerische Funktion <-- symbolische Funktion
sol = fzero(f_hand, x0)        
disp(['Loesung  ',num2str(sol),'  bei Startloesung ', num2str(x0)]);

%% Loese nach einer bestimmten Variablen symbolisch auf ("Umstellen nach")
% Stelle y = 2*x^2+4  nach x um
clear; 
%clc;
disp('Loese nach einer bestimmten Variablen auf')
syms x y;
f  = 2*x^2+4-y       %  Notwendig:  f(x,y) = 0
aa = solve(f==0,x)
pretty(aa)            % Schoenere Ausgabe

%% Schraenke das Loesen auf x>0 ein.
assume(x,'positive')
a0 = solve(f==0,x)

%% Algebraische Gleichungssysteme
%% Lineares System
%                               a +   b = 3
%                               a + 2*b = 6
%%  Variante 1: symbolisch  (mit Parameter)
%      via      f_1(a,b) := a +   b - 3 = 0
%               f_2(a,b) := a + 2*b - c = 0 
clear; 
%clc;
disp('Lineares System  Variante 1: symbolisch')
syms a b c;
%f(1) = a +   b - 3;               % nicht in octave
%f(2) = a + 2*b - c;               % nicht in octave
% [A,B] = solve(f(1),f(2),a,b)     % nicht in octave 
[A,B] = solve(a +   b - 3==0,a + 2*b - c==0,a,b)         % 2 symb. Variable als Ergebnis
%  bzw
SOLUTION = solve(a +   b - 3==0,a + 2*b - c==0,a,b)      % Cell Array als Ergebnis
SOLUTION.a
SOLUTION.b

%% Variante 2: numerisch
%     via loesen von           K*x = f     (x entspricht [a b] aus symbolischer Loesung)
clear; 
%clc;
disp('Lineares System  Variante 2: numerisch')
K = [1 1; 1 2];           % Koeffizientenmatrix
f = [ 3; 6];              % rechte Seite
x = K\f                   % Loesen des linearen Gleichungssystems

%% Nichtlineares System
%                              3a^2 + b^2 = 1
%                               a   + b   = 1
%%  Variante 1: symbolisch
%      via      f_1(a,b) := 3a^2 + b^2 - 1 = 0
%               f_2(a,b) := a   + b   - 1 = 0 
clear; 
%clc;
disp('Nichtlineares System  Variante 1: symbolisch ')
syms a b;
%f(1) = 3*a^2 + b^2 - 1;               % nicht in octave
%f(2) =   a   + b   - 1;               % nicht in octave
%[A,B] = solve(f(1),f(2))              % nicht in octave
[A,B] = solve(3*a^2 + b^2 - 1==0,a   + b   - 1==0)
%  bzw
SOLUTION = solve(3*a^2 + b^2 - 1==0,a   + b   - 1==0)
%SOLUTION.a         % nicht in octave % beide Loesungen werden ausgegeben
%SOLUTION.b         % nicht in octave
% Eine einzelne Loesung extrahieren
%Solution = [SOLUTION.a, SOLUTION.b]    % nicht in octave
Solution = SOLUTION{1};                 % matlab:  SOLUTION(1)
Solution.a
Solution.b

%% nochmal symbolisch, aber mit Parameter in Loesung
%                              3a^2 + b^2 = c
%                               a   + b   = c^2
disp('Nichtlineares System  Variante 1: symbolisch mit Variablen in Loesung')
syms c d;
%f(1) = 3*a^2 + b^2 - c;
%f(2) =   a   + b   - c^2;
%[A,B] = solve(f(1),f(2),a,b)
[A,B] = solve(3*a^2 + b^2 - c==0,a   + b   - c^2==0,a,b)

%% numerische Loesung via fsolve() [benoetigt die Optimization Toolbox];
%                              3a^2 + b^2 = 1
%                               a   + b   = 1
disp('Nichtlineares System  Variante 2: numerisch')

x_0 = [0.3, 0.6];
sol = fsolve(@(x)[3*x(1)^2 + x(2)^2 - 1; x(1)   + x(2)   - 1], x_0)
disp(['Numerische Loesung  ',num2str(sol),' bei Startloesung ',num2str(x_0)])


%% Integrieren und Differenzieren nach einer Variablen
% $$f = x^2 - 3*x + 4$$
clear; 
%clc;

disp('Integrieren und Differenzieren nach einer Variablen')
syms x;
f = x^2 - 3*x + 4
% Differenzieren
diff(f)              % oder  diff('x^2 - 3*x + 4')

% Integrieren
int(f)               % unbestimmtes Integral
int(f,x)
int(f,0,1)           % bestimmtes Integral
int(f,x,0,1)

%% Int+Diff+Vereinfachen
syms x
f  = exp(x)*cos(x)
F = int(f,x)
Fx = diff(F,x)
Fx-f                 % das soll == 0 sein?
simplify(Fx-f)       % OK, jetzt glaube ich das.

%% Integrieren und Differenzieren mehrdimensionaler Funktionen
% $$f = [a^2 + b^2 - 1, a + b - 1]$$
clear; 
%clc;

disp('Integrieren und Differenzieren mehrdimensionaler Funktionen')
syms a b
f = [a^2 + b^2 - 1, a + b - 1]

Jac = jacobian(f);

f1_a = diff(f(1),a)            % 1. Ableitung 1. Funktion bzgl. a
f1_a = diff(f(1),2)            % 2. Ableitung 1. Funktion bzgl. a

int(f(1))                      % unbestimmtes Integral bzgl. letzer Variable

int(f(1),a)                    % unbestimmtes Integral bzgl. Variable a

int(f(1),a,1,2)                % bestimmtes Integral bzgl. Variable a 

%% Gradient, Laplace, Rotation, Divergenz
%
%
clear
%clc
syms x y z a
f = x^3+sin(y)^2+exp(a*z)

% Gradient von f
disp('Gradient von f')
gradient(f,[x,y,z])

% Laplace von f
disp('Laplace von f')
laplacian(f,[x,y,z])

% Divergenz von g
g = [x^3+sin(y)^2+exp(a*z); exp(x)-log(z); y*z+x*y]
disp('Divergenz von g')
divergence(g,[x,y,z])

% Rotation von g
disp('Rotation von g')
curl(g,[x,y,z])

%% Grenzwerte:
% from docu for limit  :
% $$ \lim_{x\to\infty} (1+\frac{a}{x})^x $$
syms x a;
pretty(limit((1 + a/x)^x, x, inf))

 
%% Evaluate the sum of the following multivariable expression with respect to k:
%  from docu for 'symsum'
%
% $$\sum_{k=0}^\infty \frac{x^k}{k!}$$
%
syms x k;
pretty(symsum(x^k/sym('k!'), k, 0, inf))

% anderes Bsp:
pretty(symsum(1/k^2, k, 1, inf))

%%
% <https://de.mathworks.com/help/matlab/matlab_prog/marking-up-matlab-comments-for-publishing.html Publishing>
% publish('v_6_b.m')




