%% Övning 2
% Table of Contents
% Uppgift 2:
%     b) 2a ordn. diff.ekv. skriven som system av ode.
% 3. Matematisk pendel
% 4. Kula på gummiband rör sig i xy-planet.
% 5. Finite differences for boundary-value problems 1
% 6. Finite differences for boundary-value problems 2
%% Uppgift 2:
%% b) 2a ordn. diff.ekv. skriven som system av ode.
% Mass m on a non-linear spring with damper moves freely along one axis.
% Find displacement over time given starting displacement x0.

close all

% Problem parameters
x0 = 1;
T = 100;
k = 0.5; %spring constant
c = 0.01; % Damper constant
m = 0.1; % Mass of object
tspan = [0 T];

% The diff equation as system of ODEs.
f = @(t, U, k, m, c) [U(2); -k/m*U(1)^3 - c/m*U(2)];
u0 = [x0; 0];

% Solver
[t, u] = ode45(@(t, U)f(t, U, k, m, c), tspan, u0);

% Plot displacement 
plot(t, u(:,1))
grid on
xlabel('t')
ylabel('x [displacement]')
title('Displacement of frictionless mass on non-linear spring with damper')

%% 3. Matematisk pendel
% We approximate the pendulum with
% \sin \left(\varphi \right)\approx \varphi for small values (approx -0\ldotp 5<\varphi <0\ldotp 5).
% Thus,
% {\varphi^{\prime } }^{\prime } =-\frac{g}{L}\sin \left(\varphi \right)\to {\varphi^{\prime } }^{\prime } =-\frac{g}{L}\varphi
% This can be solved analytically (see Stewart 17.1)
%
% In brief: this is a homogenous differential 2nd order linear differential equation of the type{{\textrm{ay}}^{\prime } }^{\prime } +{\textrm{by}}^{\prime } +\textrm{cy}=0 with constant coeffs a,b,c.
%      This has the characteristic equation {\textrm{ar}}^2 +\textrm{br}+c=0. In this case the roots are complex r_1 ,r_2 =\alpha \pm i\beta and the general solution is y=e^{\alpha \;x} \left(c_1 \cos \left(\beta x\right)+c_2 \sin \left(\beta x\right)\right).
%     Here, we clearly have \alpha =0,\;\;\;\beta =\sqrt{\frac{g}{L}}.
% Thus, y=c_1 \cos \left(\beta x\right)+c_2 \sin \left(\beta x\right)\; and y^{\prime } ={-c}_1 \beta \sin \left(\beta x\right)+c_2 \beta \cos \left(\beta x\right).
% Initial values y\left(0\right)=c_1 =\varphi_0 and  y^{\prime } \left(0\right)=c_2 \beta =0.
%
% Putting it all together:
% \varphi \left(t\right)=\varphi_0 \cos \left(\sqrt{\frac{g}{L}}t\right)
% Using the relationship 2\pi f=\sqrt{\frac{g}{L}} we can obtain the period T=\frac{1}{f}=2\pi \sqrt{\frac{L}{g}}.

L = 0.1;
g = 9.81;

T = 2*pi*sqrt(L/g);

%
% This gives T=0.6344, which is approximately correct for values of \varphi <0\ldotp 5\;\textrm{rad}.
% Compare with painstakingly computed full differential solution below (from lab2, uppgift 5):
%
% The point of this exercise is to show that, for certain problem domains, the problem can be greatly simplified through approximations. Here we forego any numerical computation entirely.
%     Indeed sin(x)=x is the first order taylor series approximation. A higher order approximation could have been used, but the equation would have become non-linear.
%% 4. Kula på gummiband rör sig i xy-planet.
%

close all 
m = 0.12;
L = 0.4;
k = 10;
g = 9.81;
u0 = [0.3; 0; -0.4; 0]; % Begynnelsevärden
T = 10;
tspan = [0 T];
[t, u] = ode45(@(t,U)uprim(t,U,L,m,g,k), tspan, u0);

% Plot
scatter(u(:,1),u(:,3),25, linspace(1,10,length(t)), 'filled')
hold on
plot(u(1,1),u(1,3), 'ro') % The starting position.
axis([-0.3 0.3 -0.7 0])
title('Elastic band mass, scatterplot trajectory (initial pos highlighted red)')

plot(u(:,1),u(:,3),'b-')
axis square

%
%% 5. Finite differences for boundary-value problems 1
%

close all

g0 = 1; % Boundary conditions.
gL = 0;

L = pi/2; % Size of problem domain.
N = 30; % Number of intervals.
h = L/(N+1); % dx
i = 0:N+1; % Enumerated grid points.
x = i.*h; % x-positions
x_inner = x(2:end-1); % End points are excluded from computations since they are fixed.

% Set up matricies
ones_vec = ones(N,1);
A = spdiags([ones_vec, (h^2-2).*ones_vec, ones_vec], [-1 0 1],N,N); % See lab2 for more details.
full(A);
b = h.^2.*cos(pi.*x_inner)'; % b depends on x-position.
% Include the constant endpoint values in boundaries of "inner" grid. Remember, they are
% part of the l.h.s. expression that we have moved over to vector b.
b(1) = b(1) - g0;
b(end) = b(end) - gL;
u = A\b; % Solution
u = [g0; u; gL]; % Add end-points.

plot(x,u)
grid on
title('Finite differences, boundary-value problem 1')
xlabel('x');ylabel('u(x)')

%
%% 6. Finite differences for boundary-value problems 2
%
% OBS! In this case, the boundary conditions are different derivatives!
% (Neumann boundary, instead of Dirichlet).
% We approximate u'(0) on the boundary with u^{\prime } \left(0\right)=\frac{\left(u_1 -u_0 \right)}{h}=1.
% Thus, u_0 =u_1 -h. We substitute this into the expression for f_1: u_2 +u_1 \left(h^2 \cos^2 x_1 -1\right)=h^2 f\left(x_1 \right)+h
% Finally, when inserting the boundary values, we use the expression for u_0 above.

close all

g0 = 1; % Boundary conditions.
gL = 0;

L = 1; % Size of problem domain.
N = 10000; % Number of intervals.
h = L/(N+1); % dx
i = 0:N+1; % Enumerated grid points.
x = (i.*h)'; % x-positions. Column vector.
x_inner = x(2:end-1); % End points are excluded from computations since they are fixed.

f = @(x) x.*exp(x); % r.h.s function.

% Set up matricies
ones_vec = ones(N,1);
A = spdiags([ones_vec, (h^2*cos(x_inner).^2-2).*ones_vec, ones_vec], [-1 0 1],N,N); % See lab2 for more details.
A(1,1) = (h^2*cos(x_inner(1)).^2-1);
full(A);
b = h.^2.*f(x_inner); % b depends on x-position.
% Include the constant endpoint values in boundaries of "inner" grid. Remember, they are
% part of the l.h.s. expression that we have moved over to vector b.
b(1) = b(1) + h;
b(end) = b(end) - gL;
u = A\b; % Solution
u = [u(1)-h; u; gL]; % Add end-points.

plot(x,u)
grid on
title('Finite differences, boundary-value problem 2: Neumann bound')
xlabel('x');ylabel('u(x)')

%

function udot = uprim(t, U, L, m, g, k) % Uppgift 4
    r = sqrt(U(1).^2 + U(3).^2);
    if (r >= L)
        s = k*(r-L);
    else
        s = 0;
    end

    udot = [U(2); -s/r*U(1); U(4); -s/r*U(3)-m*g];
end