We are going to solve the differential equation with the boundary conditions
Let’s take the simpler boundary condition
.
Now, the first thing which we wanna do is to convert this equation in the first order form. After that, we can simply use the 4th order Runge-Kutta to obtain higher order solution.
To transform our 2nd order differential equation into first order, we define two variables. To transform nth order equations we need n such variables. This is the standard process to convert any differential equation into the system of first order differential equations.
Lets take
(1)
(2) [derivative of
]
So, we have transformed the 2nd order differential equation into two first order differential equations.
And our boundary condition will become
Now for the differential equations of the type
have the solution (for positive
). Now this solution doesn’t satisfy our condition
. Exponential functions are monotonic and hence cannot be 0 at two points. So, this tells us that
is not positive, i.e.,
.
Now, if , then our eqn becomes
The solution to this is of the form of cosines and sines, which has potential to satisfy the boundary condition. So, we have .
Defining Function in MATLAB:
function [ rhs ] = yfunc( x,ic,dummy,beta ) y1=ic(1); y2=ic(2); rhs = [y2; (beta -100)*y1]; end
Calculating Solutions:
clear; close all; clc set(0,'defaultlinelinewidth',2); set(0,'defaultaxeslinewidth',1.5); set(0,'defaultpatchlinewidth',2); set(0, 'DefaultAxesFontSize', 14) tol=10^(-4); % define a tolerance level to be achieved % by the shooting algorithm xspan=[-1 1]; %boundary conditions % define the span of the computational domain A=1; % define the initial slope at x=-1 ic=[0 A]; % initial conditions: x1(-1)=0, x1?(-1)=A n0=100; % define the parameter n0 dbeta=1; betasol=[]; beta_start=n0; % beginning value of beta for jj=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('yfunc',xspan, ic,[],beta); % solve ODEs if y(end,1)*((-1)^(jj+1)) > 0 %this checks if the beta needs to be higher or lower for convergence beta = beta-dbeta; else beta = beta + dbeta/2; %uses bisection to converge to the solution dbeta=dbeta/2; end if abs(y(end,1)) < tol % check for convergence betasol=[betasol beta]; % write out eigenvalue break % get out of convergence loop end end 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)); hold on end betasol legend(sprintf('\\beta_1 = %.4f',betasol(1)),sprintf('\\beta_2 = %.4f',betasol(2)),sprintf('\\beta_3 = %.4f',betasol(3)),... sprintf('\\beta_4 = %.4f',betasol(4)),sprintf('\\beta_5 = %.4f',betasol(5))); xlabel('x'); ylabel('\psi (x)','FontWeight','bold') saveas(gcf,'DiffEqnSol','jpg')