% EB_Beam_Shooting_Example.m 
%
% sample code demonstrating the use of the shooting method for solving an
%       eigenvalue problem, demonstrated on the Euler-Bernoulli beam equation
% problem: 
%       I*E/rho/A*u''''(x) = omega^2*u(x), 
% on the domain: 
%       a <= x <= b
% with boundary conditions
%       u(a) = 0, u'(a) = 0 (fixed end)
%       u''(b) = 0, u'''(b) = 0 (free end)
% 
% the boundary value problem is solved as an initial value problem with the
% introduction of new variables
% u1 := u
% u2 := u'
% u3 := u''
% u4 := u'''
% with these, the new eigenvalue problem is
% u1'(x) = u2(x)
% u2'(x) = u3(x)
% u3'(x) = u4(x)
% u4'(x) = u1(x)*omega^2*rho*A/I/E = u1(x)*omega^2/matpar;
% with initial condition: u1(a) = 0, u2(a) = 0, u3(a) = 0.1, u4(a) = unknown;
%   u3(a) = 0.1 is the normalization condition (the eigenvalue multiplied by
%   a real number is also an eigenvalue), while the last initial condition is
%   unknown, we have to solve the problem for this
% this can be interpreted as a nonlinear equation: the input is [omega,u4(a)], the
% output is [u3(b),u4(b)]; 
%
% the shooting method is similar in the case of the flow stability
% equations, except in the case of flow stability, the full solution is known at one end, and alpha
% (comlex wavenumber) us tuned so that the boundary condition is satisfied
% at the other end
%%
clear all
close all
clc

%% parameters
% domain boundaries
a = 0;
b = 1;
domain = [a,b];

% options for the ODE solver: absolute and relative error tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% material properties
ly = 1e-2; % heigth of the beam
lz = 1e-2; % transversal size of the beam
A = ly*lz; % area
I = lz*ly^3/12; % second moment of area
rho = 7800; % density
E = 200*1e9; % Elasticity modulus

materialpar = I*E/rho/A; % combining parameters into a single constant

%% function handles
% with these additional parameters can be passed when calling a function
%       the first one is a handle for an ordinary differential equation, which
%       can be passed to an ODE solver, which is a build-in MATLAB function for 
%       solving an Ordinary Differential Equation
% the second one is a nonlinear equation, which is the solution of the
%       Euler-Bernoulli equation as an initial value problem
%
% first function handle: ODE function
% the parameters of the ODE function handle are 
%       the independent variable (x) and
%       the depedent (u(x)) variables
% here we have two additional parameters, the constant containing the
%       material properties (matpar), the other is omega
% to be able to solve the problem, this function must be passed to the
%   nonlinear equation solver, therefore we create a function handle 
fhandle = @(x,u,matpar,omega) EBbeamODE(x,u,matpar,omega);
%
% second function handle: nonlinear equation to be solved
% the parameters of the nonlinear equation function are the boundary value
%       problem formulated as a function handle ODE (fhandle), 
%       the the domain of the computation (domain),
%       the parameter containing the two unknowns ([omega,u40]),
%       the normalization condition (u'''(a) = u30), which is chosen to be 0.1
%       and the material parameters (matpar)
% the nonlinear equation solver accepts a function that has only one
%       input parameter, so we have to 'burn in' the rest
%
u30 = 0.1;
ftosolve = @(solpar) EBbeamShooting(fhandle,domain,solpar,u30,materialpar);

%% initial guess for the nonlinear problem
input0 = [8*2*pi,0.1];
% input0 = [50*2*pi,0.1];
% input0 = [150*2*pi,0.1];
% input0 = [900*2*pi,0.1];
% input0 = [1700*2*pi,0.1];
% input0 = [3000*2*pi,0.1];

%% solution of the nonlinear problem - the unknown first derivative
input_sol = fsolve(ftosolve,input0);

%% recalculating the solution of the Euler-Bernoulli equation with the known eigenvalue and boundary (initial) condition
[x,sol] = ode45(@(x,u)fhandle(x,u,materialpar,input_sol(1)),[a,b],[0;0;u30;input_sol(2)],odeopt);

%% plotting solution
figure; 
plot(x,sol(:,1),'LineWidth',1.3);
hold on;
% derivatives
% plot(x,sol(:,2),'LineWidth',1.3);
% plot(x,sol(:,3),'LineWidth',1.3);
% plot(x,sol(:,4),'LineWidth',1.3);
xlabel('x');
ylabel('u(x)');
title(['\omega = ',num2str(input_sol(1)),', f = ', num2str(input_sol(1)/2/pi)]);

%% functions for the nonlinear problem and the ODE

function out = EBbeamShooting(fhandle,domain,solpar,u30,matpar)
% function formulating the nonlinear problem of solving the boundary value
%       problem as an initial value problem
% 
% parameters: 
%       fhandle: function handle for the ODE solver
%       domain: domain for ODE solution
%       solpar: [omega, u40] - variables to solver for
%       u30: the second derivative of u30
%       matpar: material parameters
% 
% in order to solve with fsolve, the additional inputs 
%       (fhandle, domain, u30, matpar) need to be added calling with fsolve
%       this can be done by creating a function handle

% ode option object for ODE tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% initial condition for the ODE
Icond = [0;0;u30;solpar(2)];

% creating a function handle for the ODE solver, which now contains omega,
%       and the new function handle can be given to the ODE solver ode45
fh = @(x,u) fhandle(x,u,matpar,solpar(1));

% solution of the ODE with solver ode45
%       parameters: function handle, solution domain, initial condition, 
%       ode option object
[~,x_] = ode45(fh,domain,Icond,odeopt);

% the output is the boundary condition at the other end: 
%       [u3(b); u4(b)] = [u''(b); u'''(b)] ( = [0; 0])
out = x_(end,[3,4]).';
end

function dudt = EBbeamODE(x,u,matpar,omega)
% ODE function for the Poisson-equation
% u'(1) = u(2);
% u'(2) = u(3);
% u'(3) = u(4);
% u'(4) = omega^2*rho*A/I/E * u(1) = omega^2/matpar*u(1)
% parameters:
% x: independent variable
% u: dependent variable - 4x1 column vector
% matpar: constant containing the material properties
% omega: eigenvalue
% in order to use with ode solver (e.g. ode45), the inputs matpar and omega
%       need to be added 
%       this can be done by creating a function handle

dudt = zeros(4,1);

dudt(1,1) = u(2);
dudt(2,1) = u(3);
dudt(3,1) = u(4);
dudt(4,1) = omega^2/matpar*u(1);

end