Scientific Computing - Quiz 4 - Problem 3-13

Problem:

To be copied

Solution:

The differential equation is so easy so we can first solve it analytically, again, using Laplace transform.

First, let $ y(x) = \psi(x - 1) $, that allow us to pose the initial condition at $ x = -1 $ to $ x = 0 $.

$ \begin{eqnarray*} y'' &=& \psi''(x - 1) \\ &=& (\beta - 100)\psi(x - 1) \\ &=& (\beta - 100)y \\ \mathcal{L}(y'') &=& \mathcal{L}((\beta - 100)y) \\ s\mathcal{L}(y') - y'(0) &=& (\beta - 100)\mathcal{L}(y) \\ s\mathcal{L}(y') - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s(s\mathcal{L}(y) - y(0)) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s(s\mathcal{L}(y) - 0) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s^2\mathcal{L}(y) - 1 &=& (\beta - 100)\mathcal{L}(y) \\ s^2\mathcal{L}(y) + (100 - \beta)\mathcal{L}(y) &=& 1 \\ \mathcal{L}(y) &=& \frac{1}{s^2 + (100 - \beta)} \\ \mathcal{L}(y) &=& \frac{1}{\sqrt{(100 - \beta)}}\frac{\sqrt{(100 - \beta)}}{s^2 + (100 - \beta)} \\ y &=& \frac{1}{\sqrt{(100 - \beta)}}\sin(\sqrt{(100 - \beta)}x) \\ \psi(x) &=& \frac{1}{\sqrt{(100 - \beta)}}\sin(\sqrt{(100 - \beta)}(x + 1)) \\ \end{eqnarray*} $

Here is the MATLAB code for solving the ODE:

clear all; close all; clc

main driver (just paste into MATLAB)
tol = 10^(-4);
xspan = [-1 1];
A = 1;
ic = [0 A];
beta = 91;
[t y] = ode45('f_rhs', xspan, ic, [], beta);
plot(t, y(:,1))
y(end, 1)

and the function to be solved:

function rhs=f_rhs(x, ic, dummy, beta)
  y1 = ic(1);
  y2 = ic(2);
  rhs = [y2;
         (beta - 100)*y1];

this is the overall code for finding all eigenvalues

clear all;                                         % clear all previously defined variables
close all;                                         % clear all previously defined figures

tol=10^(-4);                                       % define a tolerance level to be achieved by the shooting algorithm
col=[’r’,’b’,’g’,’c’,’m’,’k’];                     % eigenfunction colors
n0=100;                                            % define the parameter n0
A=1;                                               % define the initial slope at x=-1
x0=[0 A];                                          % initial conditions: x1(-1)=0, x1’(-1)=A
xp=[-1 1];                                         % define the span of the computational domain
beta_start=n0;                                     % beginning value of beta

for modes=1:5                                      % begin mode loop
  beta=beta_start;                                 % initial value of eigenvalue beta
  dbeta=n0/100;                                    % default step size in beta
  for j=1:1000                                     % begin convergence loop for beta
    [t,y]=ode45(’f_rhs’,xp,x0,[],n0,beta);         % solve ODEs
    if abs(y(end,1)-0) < tol                       % check for convergence
      beta                                         % write out eigenvalue
      break                                        % get out of convergence loop
    end
    if (-1)^(modes+1)*y(end,1)>0                   % this IF statement block
      beta=beta-dbeta;                             % checks to see if beta
    else                                           % needs to be higher or lower
      beta=beta+dbeta/2;                           % and uses bisection to
      dbeta=dbeta/2;                               % converge to the solution
    end %
  end % end convergence loop
  beta_start=beta-0.1;                             % after finding eigenvalue, pick new starting value for next mode
  norm=trapz(t,y(:,1).*y(:,1))                     % calculate the normalization
  plot(t,y(:,1)/sqrt(norm),col(modes)); hold on    % plot modes
end                                                % end mode loop

Now we can check with the analytic answers - note the output beta values actually lead to full/half cycle of the sine.

In order for $ \psi(1) = 0 $, we need $ \frac{1}{\sqrt{(100 - \beta)}}\sin(2\sqrt{(100 - \beta)}) = 0 $, of course that correspond to $ \sqrt{100 - \beta} $ is a multiple of $ \frac{\pi}{2} $