%% Solution of Assignment 2 (MATH246 Spring 2024)
%

clearvars; clf; format short g

%% Problem 1(a)
% Solve the homogeneous problem Ly=0: Plugging in y(t)=e^{rt} gives
%        Le^{rt}=(r^2+r-2)e^rt\stackrel{!}{=}0
% Hence we need r^2+r-2=0. The quadratic formula gives r_1 =1,\; r_2=-2, yielding the fundamental set of solutions
%         Y_1(t)=e^{t},\qquad Y_2(t)=e^{-2t}
% Find a particular solution: For Y(t)=A\cos t+B\sin t we obtain
%        LY = A(-\cos t-\sin t-2\cos t)+ B(-\sin t+\cos t-2\sin t) = (-3A+B)\cos t+(-A-3B)\sin t
% This gives the linear system \left[
% \matrix{-3 & 1 \cr -1 &-3}
% \right]
% \left[
% \matrix{A\cr B}
% \right]
% =\left[
% \matrix{0\cr 1}
% \right] which has the solution

[-3 1;-1 -3]\sym([0;1])      % A\b solves linear system  A x = b

% We obtain the particular solution Y(t)=-\frac{1}{10}\cos t-\frac3{10}\sin t.
% Hence the general solution of Ly=f is y(t)=Y(t)+c_1 Y_1(t)+c_2 Y_2(t):
%         y(t) = -\frac{1}{10}\cos t-\frac3{10}\sin t+c_1 e^{t}+c_2 e^{-2t},\qquad c_1,c_2\in\mathbb{R}
% We determine c_1,c_2 from the initial conditions:
%         \left[
% \matrix{ Y_1(t_0) & Y_2(t_0) \cr
% Y_1'(t_0) & Y_2'(t_0) }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ y_0 -Y(t_0) \cr y_0'-Y'(t_0) }
% \right]  gives  \left[
% \matrix{ 1 & 1 \cr
% 1 & -2 }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ 1+\frac1{10} \cr 2 + \frac3{10}}
% \right]
% which has the following solution \left[
% \matrix{c_1 \cr c_2}
% \right]:

[1 1; 1 -2]\sym([1.1; 2.3])

% Hence  the solution of the initial value problem is
%         y(t) = -\frac{1}{10}\cos t-\frac3{10}\sin t+\frac32 e^{t}-\frac25 e^{-2t}
%% Problem 1(b)
% Note that the coefficient functions are not continuous at t=0. Since the initial condition is given at t_0=1 we will obtain solutions which exist on the interval (0,\infty), i.e., for t>0.
% Solve the homogeneous problem Ly=0: Plugging in y(t)=t^r gives
%        Lt^r=(r(r-1)+r-1)t^{r-2}
% =(r^2-1)t^{r-2}\stackrel{!}{=}0
% Hence we need r^2-1=0. We get the solutions r_1 =1,\; r_2=-1, yielding the fundamental set of solutions
%         Y_1(t)=t,\qquad Y_2(t)=t^{-1}
% Find a particular solution: For Y(t)=At^s we obtain
%        LY =A(s^2-1)t^{s-2} \stackrel{!}{=}t^2
% which gives s=4 and 15A=1. Hence we get the particular solution Y(t)=\frac1{15}t^4.
% Hence the general solution of Ly=f is y(t)=Y(t)+c_1 Y_1(t)+c_2 Y_2(t):
%         y(t) =
% \frac1{15}t^4+c_1t+c_2 t^{-1}
% ,\qquad c_1,c_2\in\mathbb{R}
% We determine c_1,c_2 from the initial conditions:
%         \left[
% \matrix{ Y_1(t_0) & Y_2(t_0) \cr
% Y_1'(t_0) & Y_2'(t_0) }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ y_0 -Y(t_0) \cr y_0'-Y'(t_0) }
% \right]  gives   \left[
% \matrix{ 1 & 1 \cr
% 1 & -1 }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ 1-\frac1{15}\cr 1-\frac4{15}}
% \right]
% which has the following solution \left[
% \matrix{c_1 \cr c_2}
% \right]:

[1 1; 1 -1]\[ 1-1/sym(15) ; 1-4/sym(15) ]

% Hence  the solution of the initial value problem is
%         y(t) =
% \frac1{15}t^4+\frac56t+\frac1{10} t^{-1}
%
%% Problem 2(a)
% We already solved the homogeneous ODE in 1(a). For the particular solution we use Y(t)=A_0t+A_1 and get
%        LY = 0+A_0 -2(A_0 t+A_1) = -2A_0 t+(A_0-2A_1)\stackrel{!}{=} 1\cdot t + 0
% This gives -2A_0=1\iff A_0=-\frac12 and A_0-2A_1=0, hence A_1=-\frac14 yielding Y(t)=-\frac{t}2-\frac14.
% We determine c_1,c_2 from the initial conditions:
%         \left[
% \matrix{ Y_1(t_0) & Y_2(t_0) \cr
% Y_1'(t_0) & Y_2'(t_0) }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ y_0 -Y(t_0) \cr y_0'-Y'(t_0) }
% \right]  gives  \left[
% \matrix{ 1 & 1 \cr
% 1 & -2 }
% \right]
% \left[
%  \matrix{c_1 \cr c_2 }
% \right]
% = \left[
%   \matrix{ 1+\frac1{4} \cr 2 + \frac1{2}}
% \right]
% which has the following solution \left[
% \matrix{c_1 \cr c_2}
% \right]:

[1 1; 1 -2]\[5/sym(4); 3/sym(2)]

% Hence  the solution of the initial value problem is
%         y(t) =-\frac{t}2-\frac14+\frac43e^t-\frac1{12}e^{-2t}
%% Problem 2(b)
% We rewrite the problem as a first order system of ODEs \vec y'=\vec f(t,\vec y):
%            y_1(t)=y(t),\quad
% y_2(t)=y'(t),\qquad
% \left[
% \matrix{y_1'(t) \cr y_2'(t)}
% \right]
% =\left[
% \matrix{y_2(t) \cr t-y_2(t)+2y_1(t)}
% \right]

f = @(t,y) [ y(2) ; t-y(2)+2*y(1) ];

% with initial condition \vec y(0) = \left[
% \matrix{1 \cr 1}
% \right]:

t0 = 0; y0 = [1;1];

% Perform 1 step of improved Euler with h=\frac12:

h = 1/2;

s1 = f(t0,y0)
yE = y0 + h*s1
s2 = f(t0+h,yE)
y1 = y0 + h*(s1+s2)/2

% We obtain for \vec y(\frac12) the approximation \left[
% \matrix{1.625 \cr 1.75}
% \right], hence y(\frac12) is approximated by 1.625.
%% Problem 2(c)
% The code for the function IEulermethod is at the end of this document.
% We want to see how the error improves when we double n, and print out
%             \text{factor}=\frac{\text{error}_{n/2}}{\text{error}_{n}}

yex = @(t) -t/2-1/4+4/3*exp(t)-1/12*exp(-2*t);   % solution from (a)
T = 2;
exact = yex(T);
olderr = NaN;
for n=2.^(3:8)
  [ts,ys] = IEulerMethod(f,[0,T],y0,n);
  err = ys(end,1)-exact;
  fprintf('n=%3g:  error=%12e,  factor=%g\n',n,err,olderr/err)
  olderr = err;
end

% The Improved Euler Method has convergence order 2: this means that the error at the final time T satisfies
%          \text{error}\le \frac{C}{n^2}
% This upper bound gets divided by 4 when we double the number n. We see that the actual errors decrease in the same way, as the computed factors get closer and closer to 4.
%% Problem 2(d)

format long g    % show 15 digits of results
T = 8;
exact = yex(T);
[ts,ys] = ode45(f,[0,T],y0);
Y = ys(end,1)
err = Y - exact
opt = odeset('RelTol',3e-14,'AbsTol',1e-14);
[ts,ys] = ode45(f,[0,T],y0,opt);
Y = ys(end,1)
err = Y - exact

%
%% Problem 3(i)(a)
% We rewrite the problem as a first order system of ODEs \vec y'=\vec f(t,\vec y):
%            y_1(t)=y(t),\quad
% y_2(t)=y'(t),\qquad
% \left[
% \matrix{y_1'(t) \cr y_2'(t)}
% \right]
% =\left[
% \matrix{y_2(t) \cr -y_2(t)+2y_1(t)}
% \right]

f = @(t,y) [ y(2) ; -y(2)+2*y(1) ];
vectorfield(f,-2:.25:2,-2:.25:2); hold on
for a=-2:2
  for b=-2:2
    for T=[8,-8]
      [ts,ys] = ode45(f,[0 T],[a;b]);  % find solution for t from 0 to T
      plot(ys(:,1),ys(:,2),'b');       % plot trajectory in the phas plane
    end
  end
end
hold off; axis([-2 2 -2 2]); xlabel('y_1'); ylabel('y_2')
title('25 trajectories for t\in [-8,8]')

%% Problem 3(i)(b)
% We found in 1(a): the general solution of the homogeneous ODE is
%           y(t) = c_1e^t +c_2 e^{-2t},    y'(t) = c_1e^t +c_2 (-2)e^{-2t}
% This gives for \vec y(t)= \left[
% \matrix{y(t) \cr y'(t) }
% \right]
%           \vec y(t) = c_1e^t \left[
%  \matrix{1 \cr 1}
% \right]
% +c_2e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
% (1) For c_1=1,\; c_2=0 we get \vec y(t) = e^t \left[
%  \matrix{1 \cr 1}
% \right]
%  for t\in (-\infty,\infty). Hence e^t goes from 0 to \infty, so the solution starts at the origin and moves along a half line going in the direction \left[\matrix{1\cr 1}\right] (dashed blue line below).
% (2) For c_1=-1,\; c_2=0 we get \vec y(t) = -e^t \left[
%  \matrix{1 \cr 1}
% \right]
% : a half line starting at the origin and going in the direction \left[\matrix{-1\cr -1}\right] (dotted blue line below).
% (3) For c_1=0,\; c_2=1 we get \vec y(t) = e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
%  for t\in (-\infty,\infty). Hence e^{-2t} goes from \infty to 0, so the trajectory is a half line going in the direction \left[\matrix{1\cr -2}\right], the solution moves along this half line towards the origin (dashed green line below)
% (4) For c_1=0,\; c_2=-1 we get \vec y(t) = -e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
%  for t\in (-\infty,\infty) (dotted green line below)

plot([0 2],[0 2],'b--'); hold on
plot([0 -2],[0 -2],'b:'); 
plot([0 2],[0 -4],'g--'); 
plot([0 -2],[0 4],'g:'); 
vectorfield(f,-2:.5:2,-2:.5:2); 
legend('(1)','(2)','(3)','(4)','Location','eastoutside')
axis([-2 2 -2 2]); hold off;
title('Four trajectories for t\in(-\infty,\infty)')

%
%% Problem 3(ii)(a)
% We rewrite the problem as a first order system of ODEs \vec y'=\vec f(t,\vec y):
%            y_1(t)=y(t),\quad
% y_2(t)=y'(t),\qquad
% \left[
% \matrix{y_1'(t) \cr y_2'(t)}
% \right]
% =\left[
% \matrix{y_2(t) \cr -\frac12y_1(t)-\frac94y_2(t)}
% \right]

f = @(t,y) [ y(2) ; -1/2*y(1)-9/4*y(2) ];
vectorfield(f,-2:.25:2,-2:.25:2); hold on
for a=-2:1:2
  for b=-2:1:2
    for T=[8,-8]
      [ts,ys] = ode45(f,[0 T],[a;b]); 
      plot(ys(:,1),ys(:,2),'b')
    end
  end
end
hold off; axis([-2 2 -2 2]); xlabel('y_1'); ylabel('y_2')
title('25 trajectories for t\in [-8,8]')

%% Problem 3(ii)(b)
% We get the characteristic equation r^2+\frac94r+\frac12=0. This gives r_1=-2,\;r_2=-\frac14, so  the general solution of the homogeneous ODE is
%           y(t) = c_1e^{-2t} +c_2 e^{-t/4},    y'(t) = c_1(-2)e^{-2t} +c_2 (-\frac14)e^{-t/4}
% This gives for \vec y(t)= \left[
% \matrix{y(t) \cr y'(t) }
% \right]
%           \vec y(t) = c_1e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
% +c_2e^{-t/4} \left[
%  \matrix{1 \cr -\frac14}
% \right]
% (1) For c_1=1,\; c_2=0 we get \vec y(t) = e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
%  for t\in (-\infty,\infty). Hence e^{-2t} goes from \infty to 0, this is a half line going in the direction \left[\matrix{1\cr -2}\right] (dashed blue line below), the solution is moving toward the origin.
% (2) For c_1=-1,\; c_2=0 we get \vec y(t) = -e^{-2t} \left[
%  \matrix{1 \cr -2}
% \right]
% : a half line starting at the origin and going in the direction \left[\matrix{1\cr -2}\right] (dotted blue line below).
% (3) For c_1=0,\; c_2=1 we get \vec y(t) = e^{-t/4} \left[
%  \matrix{1 \cr -\frac14}
% \right]
%  for t\in (-\infty,\infty). Hence e^{-t/4} goes from \infty to 0, so the trajectory is a half line going in the direction \left[\matrix{1\cr -\frac14}\right], the solution moves along this half line towards the origin (dashed green line below)
% (4) For c_1=0,\; c_2=-1 we get \vec y(t) = -e^{-t/4} \left[
%  \matrix{1 \cr -\frac14}
% \right]
%  for t\in (-\infty,\infty) (dotted green line below)

plot([0 2],[0 -4],'b--'); hold on
plot([0 -2],[0 4],'b:'); 
plot([0 8],[0 -2],'g--'); 
plot([0 -8],[0 2],'g:'); 
vectorfield(f,-2:.25:2,-2:.25:2); 
legend('(1)','(2)','(3)','(4)','Location','eastoutside')
axis([-2 2 -2 2]); hold off; 
title('Four trajectories for t\in(-\infty,\infty)')

%
%% Function IEulerMethod for problem 2

function [ts,ys] = IEulerMethod(f,interval,y,n)
% [ts,ys] = EulerMethod(f,[t0,T],y0,n)
% initial value problem y'=f(t,y), y(t0) = y0
% use n steps of Euler method to approximate solution in interval [t0,T]
t = interval(1); T = interval(2);
h = (T-t)/n;
m = size(y,1);                          % size of vector y
ts = zeros(n+1,1); ys = zeros(n+1,m);   % initialize vector for t-values, array for y-values
ts(1) = t; ys(1,:) = y.';               % y.' is transpose, gives row vector
for k=1:n
  s1 = f(t,y);
  yE = y + h*s1;
  s2 = f(t+h,yE);
  t = t + h;
  y = y + h*(s1+s2)/2;
  ts(k+1) = t; ys(k+1,:) = y.';         % store t,y in vector ts, array ys
end
end