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} $
