%% Solving differential equations using Galerkin
% Example
% \left\lbrace \begin{array}{ll}
% -\frac{d^2 u}{d\;x}=x^2  & \mathrm{for}\;x\in \left\lbrack 0,1\right\rbrack \\
% u\left(0\right)=-1 & \;\\
% \frac{d\;u}{d\;x}\left(1\right)=-2 & \;
% \end{array}\right.
% Weak form:
% \int_0^1 -\frac{d^2 u}{d\;x}\;v\;d\;x=\int_0^1 x^2 v
% \Rightarrow \int_0^1 \frac{d\;u}{d\;x}\;\frac{d\;v}{d\;x}\;d\;x=\int_0^1 x^2 v+{\left\lbrack \frac{d\;u}{d\;x}\left(x\right)\;v\left(x\right)\right\rbrack }_0^1 \;\;\;\;\;\left(1\right)
% Approximation using Galerkin
%
% \mathrm{du1}=\frac{d\;u}{d\;x}\left(0\right),\;\;\mathrm{du2}=\frac{d\;u}{d\;x}\left(1\right)=-2,\;u_1 =-1
% Approximate using one quadratic element on the whole domain
% Quadratic basis functions: For x\in \left\lbrack 0,1\right\rbrack \; we have
% \begin{array}{l}
% \varphi_1 =2\left(x-1\right)\left(x-1/2\right)\\
% \varphi_2 =-4x\left(x-1\right)\\
% \varphi_3 =2x\left(x-1/2\right)
% \end{array}
% The solution is then approximated by
% U\left(x\right)=\sum_{i=1}^3 u_i {\;\varphi }_i
% and the derivate is given by
% \frac{d\;U}{d\;x}=\frac{d}{d\;x}\left(\sum_{\;} u_{i\;} \varphi_i \left(x\right)\right)=\sum_{\;} u_i \;\frac{d{\;\varphi }_i }{d\;x}
% Insert the approximation into the weak form, replacing v\;with \varphi and \frac{d\;v}{d\;x} with \frac{d\;\varphi }{d\;x}.
% \int_0^1 \sum_{\;} \frac{d\;\varphi_i }{d\;x}\;u_i \;\frac{d\;\varphi_j }{d\;x}\;d\;x=\int_0^1 x^2 \;\varphi_j \;d\;x+\underset{-2}{\underbrace{\frac{d\;u}{d\;x}\left(1\right)\varphi \left(1\right)} }
% Switch to Matlab to solve the equations.

clear
syms u1 u2 u3 x 
phi(x) = [2*(x-1)*(x-1/2)  -4*x*(x-1)   2*x*(x-1/2)];
u = [u1; u2; u3];
dphi = diff(phi,x)


syms du0
eqs = int(dphi*u*dphi,x,0,1) == int(x^2*phi, x, 0, 1) + phi(1)*-2 - phi(0)*du0;
eqs.'

u1 = -1

% Here we have two alternatives, either reduce the equations to two, by eliminating the equation where u_1 is known, which is row 1, or by solving for \frac{d\;u}{d\;x}\left(0\right) which we denote by du0. The result is the same.

% [u2, u3, du0] = solve(subs(eqs),[u2, u3, du0] )
[u2, u3] = solve(subs(eqs(2:3)),[u2, u3] )

U = simplify(subs(phi*u))

% dU = diff(U)
% double(dU(1))

syms ue(x)
DE = -diff(ue,x,2) == x^2;
du = diff(ue,x,1);
BC = [ue(0)==-1, du(1)==-2];
ue = dsolve(DE,BC)

figure
fplot(ue, [0,1],'r','DisplayName','Exact solution'); hold on
fplot(U,[0,1],'b','DisplayName','Approximation')
plot([0,0.5,1],[u1,u2,u3],'ob','DisplayName','u_i')
legend show

%% Piecewise approximation
%
% Instead of approximating the solution using one trial function, we can split the domain into N elements and approximate the solution by finding more nodal solutions.
% In this example we want to find all 6 nodal values u_i such that
% \int_a^b \left(U-u\right)\varphi_i \;d\;x=0
%% Deriving basis functions
% In order to solve this we need to define \varphi_i \left(x\right) for each element. Let's see how basis functions are defined.
%
% We can define a polynomial function using its nodal values:
% U\left(x\right)=\frac{x_2 -x}{x_2 -x_1 }u_1 +\frac{x-x_1 }{x_2 -x_1 }=\varphi_1 u_1 +\varphi_1 u_1
% The definition of a nodal basis function: Let the i'th basis function be equal to 1 in the i'th node and 0 in all other nodes.
% \varphi_j \left(x_i \right)=\left\lbrace \begin{array}{ll}
% 1 & \mathrm{if}\;i=j\\
% 0 & \mathrm{if}\;i\not= j
% \end{array}\right.
%
% To derive the basis functions use the definition and make the ansatz:
% \varphi \left(x\right)=c_1 +c_2 x
% \begin{array}{l}
% \varphi_1 \left(0\right)=c_1 +c_2 \cdot 0=1\\
% \varphi_1 \left(1\right)=c_1 +c_2 \cdot 1=0\\
% \varphi_2 \left(0\right)=c_1 +c_2 \cdot 0=0\\
% \varphi_2 \left(1\right)=c_1 +c_2 \cdot 1=1
% \end{array}
% or on matrix form
% \left\lbrack \begin{array}{cc}
% 1 & x_1 \\
% 1 & x_2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% c_{11} \\
% c_{21}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% 1\\
% 0
% \end{array}\right\rbrack and \left\lbrack \begin{array}{cc}
% 1 & x_1 \\
% 1 & x_2
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% c_{12} \\
% c_{22}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% 0\\
% 1
% \end{array}\right\rbrack
% or
% \underset{\mathit{\mathbf{A}}}{\underbrace{\left\lbrack \begin{array}{cc}
% 1 & x_1 \\
% 1 & x_2
% \end{array}\right\rbrack } } \underset{\mathit{\mathbf{C}}}{\underbrace{\left\lbrack \begin{array}{cc}
% c_{11}  & c_{12} \\
% c_{21}  & c_{22}
% \end{array}\right\rbrack } } =\underset{\mathit{\mathbf{I}}}{\underbrace{\left\lbrack \begin{array}{cc}
% 1 & 0\\
% 0 & 1
% \end{array}\right\rbrack } } \Rightarrow \mathit{\mathbf{C}}={\mathit{\mathbf{A}}}^{-1} \mathit{\mathbf{I}}
% \varphi \left(x\right)=\left\lbrack 1\;x\right\rbrack \mathit{\mathbf{C}}

syms x1 x2 x
A = [1 x1
     1 x2]
C = A\eye(2)

phi = simplify([1 x]*C)

% For quadratic make the ansatz: \varphi \left(x\right)=c_1 +c_2 x+c_3 x^2
% \begin{array}{l}
% \varphi_1 \left(0\right)=c_1 +c_2 \cdot 0+c_3 \cdot 0^2 =1\\
% \varphi_1 \left(1/2\right)=c_1 +c_2 \cdot 1/2+c_3 \cdot {\left(1/2\right)}^2 =0\\
% \varphi_1 \left(1\right)=c_1 +c_2 \cdot 1+c_3 \cdot 1^2 =0\\
% \vdots
% \end{array}
% which leads to the matrix system
% \left\lbrack \begin{array}{ccc}
% 1 & x_1  & x_1^2 \\
% 1 & x_2  & x_2^2 \\
% 1 & x_3  & x_{\;3}^3
% \end{array}\right\rbrack \left\lbrack \begin{array}{ccc}
% c_{11}  & c_{12}  & c_{13} \\
% c_{21}  & c_{22}  & c_{23} \\
% c_{31}  & c_{32}  & c_{33}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{ccc}
% 1 & 0 & 0\\
% 0 & 1 & 0\\
% 0 & 0 & 1
% \end{array}\right\rbrack
% \varphi \left(x\right)=\left\lbrack 1\;\;x\;\;x^2 \right\rbrack \mathit{\mathbf{C}}

syms x1 x2 x3 x
A = [1 x1 x1^2
     1 x2 x2^2
     1 x3 x3^2];
C = A\eye(3);

phi = simplify([1 x x^2]*C).'

% and for x_1 =0,x_2 =1/2,x_3 =1 we have

subs(phi,[x1 x2 x3], [0, 1/2, 1])

%
%% Piecewise approximation of a differential equation using a systematic approach
% Solve
% \left\lbrace \begin{array}{ll}
% -\frac{d\;u}{d\;x}=1 & \mathrm{for}\;x\in \left\lbrack 1\;3\right\rbrack \\
% u\left(1\right)=2 & \\
% \frac{d\;u}{d\;x}\left(3\right)=-1 &
% \end{array}\right.
% using 4 linear elements.
% Solution: Weak form:
% \int_1^3 \frac{d\;u}{d\;x}\frac{d\;v}{d\;x}\;d\;x=\int_{1\;}^3 1\cdot v\;d\;x+{\left\lbrack \frac{d\;u}{d\;x}\;v\right\rbrack }_1^3
% Galerkin approximation using linear polynomials:
% U\left(x\right)=\sum_{i=1}^2 \varphi_i \left(x\right)u_i
% \varphi \left(x\right)=\left\lbrack \frac{x_2 -x}{x_2 -x_1 }\;\;\;\frac{x-x_1 }{x_2 -x_1 }\right\rbrack
% \frac{d\;u}{d\;x}\approx \frac{d\;U}{d\;x}=\frac{d}{\;d\;x}\left(\sum \varphi_i \left(x\right)u_i \right)=\frac{d\;\varphi }{d\;x}u_i
% Insert the approximation into the weak form
% \int_1^3 \sum_i \frac{d\;\varphi_i }{d\;x}u_i \;\frac{d\;\varphi_j }{d\;x}\;d\;x=\int_1^3 1\cdot \varphi_j \;d\;x+{\left\lbrack \frac{d\;u}{d\;x}\varphi_j \right\rbrack }_1^3
% =\sum_i \int_1^3 \frac{d\;\varphi_i }{d\;x}\;\frac{d\;\varphi_j }{d\;x}\;d\;x\;\;u_i =\int_1^3 \varphi_j \;d\;x+\underset{g_e }{\underbrace{\frac{d\;u}{d\;x}{\left(3\right)\varphi }_j \left(3\right)} } \underset{g_s }{\underbrace{-\frac{d\;u}{d\;x}\left(1\right)\varphi_j \left(1\right)} }
% \frac{d\;\varphi_i }{d\;x}\;\frac{d\;\varphi_j }{d\;x} leads to \left\lbrack \begin{array}{c}
% \frac{d\;\varphi_1 }{d\;x}\\
% \frac{d\;\varphi_2 }{d\;x}
% \end{array}\right\rbrack \left\lbrack \begin{array}{cc}
% \frac{d\;\varphi_1 }{d\;x} & \frac{d\;\varphi_2 }{d\;x}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% \varphi_1^{\prime } \varphi_1^{\prime }  & \varphi_1^{\prime } \varphi_2^{\prime } \\
% \varphi_2^{\prime } \varphi_1^{\prime }  & \varphi_2^{\prime } \varphi_2^{\prime }
% \end{array}\right\rbrack
% S_{i\;j} \;u_i =f_j +g_j
% 4 linear elements will lead to 5 equations to be solved for u_i ,i=1,2,\ldotp \ldotp \ldotp ,5
%
% An element e between nodes x_i and x_j contributes only to equations i\;and j.
% Example: Element stiffness matrix for element 3 between node 3 and 4:
%
% \begin{array}{l}
% \;\;\;\;\begin{array}{ccccc}
% \;1 & \;2 & \;\;3 & \;\;\;4 & \;\;5
% \end{array}\\
% \begin{array}{c}
% 1\\
% 2\\
% 3\\
% 4\\
% 5
% \end{array}\underset{\mathit{\mathbf{S}}}{\underbrace{\left\lbrack \begin{array}{ccccc}
%  &  &  &  & \\
%  &  &  &  & \\
%  &  & S_{33}  & S_{34}  & \\
%  &  & S_{43}  & S_{44}  & \\
%  &  &  &  &
% \end{array}\right\rbrack } }
% \end{array}
% \begin{array}{l}
% S_{33} =\int_{x_3 }^{x_4 } \frac{d\;\varphi_3 }{d\;x}\;\frac{d\;\varphi_3 }{d\;x}\;d\;x\\
% S_{34} =\int_{x_3 }^{x_4 } \frac{d\;\varphi_3 }{d\;x}\;\frac{d\;\varphi_4 }{d\;x}\;d\;x\\
% S_{43} =\int_{x_3 }^{x_4 } \frac{d\;\varphi_4 }{d\;x}\;\frac{d\;\varphi_3 }{d\;x}\;d\;x\\
% S_{44} =\int_{x_3 }^{x_4 } \frac{d\;\varphi_4 }{d\;x}\;\frac{d\;\varphi_4 }{d\;x}\;d\;x
% \end{array}
% We can use local node numbering to simplify things:
%
% Now, for each element e we define x_1 and x_2, compute \varphi and \frac{d\;\varphi }{d\;x} and then compute {\mathit{\mathbf{S}}}^e and {\mathit{\mathbf{f}}}^e
% Element stiffness matrix
%  {\mathit{\mathbf{S}}}^e =\int_{x_1 }^{x_2 } \frac{d\;\varphi^T }{d\;x}\frac{d\;\varphi }{d\;x}\;d\;x
% or
% {\mathit{\mathbf{S}}}^e =\left\lbrack \begin{array}{cc}
% \int_{x_1 }^{x_2 } \frac{d\;\varphi_1 }{d\;x}\frac{d\;\varphi_1 }{d\;x}\;d\;x & \int_{x_1 }^{x_2 } \frac{d\;\varphi_1 }{d\;x}\frac{d\;\varphi_2 }{d\;x}\;d\;x\\
% \int_{x_1 }^{x_2 } \frac{d\;\varphi_2 }{d\;x}\frac{d\;\varphi_1 }{d\;x}\;d\;x & \int_{x_1 }^{x_2 } \frac{d\;\varphi_2 }{d\;x}\frac{d\;\varphi_2 }{d\;x}\;d\;x
% \end{array}\right\rbrack
% and element load vector
% {\mathit{\mathbf{f}}}^e =\int_{x_1 }^{x_2 } \varphi^T f\;d\;x
% or
% {\mathit{\mathbf{f}}}^e =\left\lbrack \begin{array}{c}
% \int_{x_1 }^{x_2 } f\varphi_1 \;\mathrm{dx}\\
% \int_{x_1 }^{x_2 } f\varphi_2 \;\mathrm{dx}
% \end{array}\right\rbrack
%% Assembly process
% For each element e, we add {\mathit{\mathbf{S}}}^e and {\mathit{\mathbf{f}}}^e to the correct place in \mathit{\mathbf{S}} and \mathit{\mathbf{f}}.
%
% element 1
% element 2
% element 3
% element 4
% \begin{array}{l}
% S_{\textrm{ii}}^{\textrm{new}} =S_{\textrm{ii}}^{\textrm{old}} +S_{11}^e \\
% S_{\textrm{ij}}^{\textrm{new}} =S_{\textrm{ij}}^{\textrm{old}} +S_{12}^e \\
% \textrm{etc}
% \end{array}
% \begin{array}{l}
% f_i^{\textrm{new}} =f_i^{\textrm{old}} +f_1^e \\
% f_j^{\textrm{new}} =f_j^{\textrm{old}} +f_2^e
% \end{array}
%
%
% Contribution from element 1
% {\mathit{\mathbf{S}}\left(\left\lbrack 1,2\right\rbrack ,\left\lbrack 1,2\right\rbrack \right)}_{2\times 2} ={\mathit{\mathbf{S}}\left(\left\lbrack 1,2\right\rbrack ,\left\lbrack 1,2\right\rbrack \right)}_{2\times 2} +S_{2\times 2}^{e_1 }
% \begin{array}{l}
% \;\;\;\;\begin{array}{ccccc}
% \;1 & \;\;\;2 & \;\;\;3 & \;\;4 & 5
% \end{array}\\
% \begin{array}{c}
% 1\\
% 2\\
% 3\\
% 4\\
% 5
% \end{array}\left\lbrack \begin{array}{ccccc}
% S_{11}^1  & S_{12}^1  &  &  & \\
% S_{21}^1  & S_{22}^1  &  &  & \\
%  &  &  &  & \\
%  &  &  &  & \\
%  &  &  &  &
% \end{array}\right\rbrack
% \end{array}
% Contribution from element 2
% {\mathit{\mathbf{S}}\left(\left\lbrack 2,3\right\rbrack ,\left\lbrack 2,3\right\rbrack \right)}_{2\times 2} ={\mathit{\mathbf{S}}\left(\left\lbrack 2,3\right\rbrack ,\left\lbrack 2,3\right\rbrack \right)}_{2\times 2} +S_{\;}^{e_2 }
% \begin{array}{l}
% \;\;\;\;\begin{array}{ccccc}
% \;1 & \;\;\;\;\;\;\;\;\;\;\;2 & \;\;\;\;\;\;\;\;\;\;3 & \;\;4 & 5
% \end{array}\\
% \begin{array}{c}
% 1\\
% 2\\
% 3\\
% 4\\
% 5
% \end{array}\left\lbrack \begin{array}{ccccc}
% S_{11}^1  & S_{12}^1  &  &  & \\
% S_{21}^1  & S_{22}^1 +S_{11}^2  & S_{12}^2  &  & \\
%  & S_{21}^2  & S_{22}^2  &  & \\
%  &  &  &  & \\
%  &  &  &  &
% \end{array}\right\rbrack
% \end{array}
% etc
%% Boundary conditions
%
%
% Since u_1 is known, we remove the associated equation from the system.
% We also need to modify the right-hand side, since u_1 \;is non-zero it will cause a reaction force that need to be added to the force vector.
% {\mathit{\mathbf{f}}}^{\mathrm{new}} ={\mathit{\mathbf{f}}}^{\mathrm{old}} +\mathit{\mathbf{g}}-\left\lbrack \begin{array}{c}
% S_{11} \\
% S_{21} \\
% \vdots \\
% S_{51}
% \end{array}\right\rbrack u_1
% and finally our equation system is
% \left\lbrack \begin{array}{cccc}
% S_{22}  & S_{23}  & S_{24}  & S_{25} \\
% S_{32}  &  &  & \vdots \\
% S_{42}  &  &  & \vdots \\
% S_{52}  & \cdots  & \cdots  & S_{55}
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% u_2 \\
% u_3 \\
% u_4 \\
% u_5
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% f_2 +0-S_{21} u_1 \\
% f_3 +0-S_{31} u_1 \\
% f_4 +0-S_{41} u_1 \\
% f_5 -2-S_{51} u_1
% \end{array}\right\rbrack
%
%
%% Matlab implementation
% \left\lbrace \begin{array}{ll}
% -\frac{d\;u}{d\;x}=1 & \mathrm{for}\;x\in \left\lbrack 1\;3\right\rbrack \\
% u\left(1\right)=2 & \\
% \frac{d\;u}{d\;x}\left(3\right)=-1 &
% \end{array}\right.

clear
a = 1; b = 3;
N = 4;
xnod = linspace(a,b,N+1);
nodes = [(1:N)',(2:N+1)'];

syms x
Phi = @(x1,x2,x)[(x-x2)/(x1-x2), -(x-x1)/(x1-x2)];
Dphi = matlabFunction(diff(Phi,x));

neq = N+1;
S = zeros(neq,neq); F = zeros(neq,1);

for iel = 1:N
    inods = nodes(iel,:);
    x1 = xnod(inods(1));
    x2 = xnod(inods(2));

    phi = Phi(x1,x2,x);
    dphi = Dphi(x1,x2);

    Se = double(int(dphi.'*dphi,x,x1,x2));
    fe = double(int(1*phi.',x,x1,x2));
    
    S(inods,inods) = S(inods,inods) + Se;
    F(inods) = F(inods) + fe;
end

presc = 1;
alldofs = 1:N+1;
free = setdiff(alldofs,presc);

u = zeros(neq,1); u(presc) = 2;
g = zeros(neq,1); 
g(1) = -0;
g(end) = -1;

F = F+g-S(:,presc)*u(presc);

u(free) = S(free,free)\F(free)


syms ue(x)
DE = -diff(ue,x,2) == 1;
du = diff(ue,x,1);
BC = [ue(1)==2, du(3)==-1];
ue = dsolve(DE,BC)

figure
fplot(ue, [1,3],'r','DisplayName','Exact solution'); hold on
plot(xnod,u,'ob-','DisplayName','Approximation')
legend show

% Computing the error
% L_2 \mathrm{err}:=\sqrt{\int_1^3 {\left(U-u_e \right)}^2 \;d\;x}

L2err = 0;
for iel = 1:N
    inods = nodes(iel,:);
    x1 = xnod(inods(1));
    x2 = xnod(inods(2));
    phi = Phi(x1,x2,x);

    U = phi*u(inods);

    L2err = L2err + int((U-ue)^2,x,x1,x2);
end
L2err = double(sqrt(L2err))

%
%