% Ex2abstractrootfinding
% Solve a 2 point Boundary Value Problem by the shooting method
% Cast the problem y'' + 4y' - 6y = x, y(0)=1, y(2)=6 as an
% initial value problem:
% Let z1 = y and z2 = y'. The problem is converted to a system:
% z1' = z2
% z2' = x - 4z2 + 6z1
% z1(0) = 1 but z2(0) is unknown. Find the correct value,
% called A here, by setting up a function and finding its root
% The function is
% rightBC(A) = [y(x=2) when y'(0)=A] - 6
% We have to solve the system to find [y(x=2) when y'(0)=A] so
% this is not a standard function definition
% to allow use to nest functions, make this script a function file:
function Ex2abstractrootfinding
close all
clear all
hold off
% set grid parameters for the domain:
dx = 1e-3;
x = 0:dx:2;
% find a bracketing interval to minimize the function rightBC
% one A value where the function is negative, and one where the
% function is positive.
Amin = -10;
t = rightBC1(Amin);
fprintf('%e: function when Amin=%f\n',t,Amin)
Amax = 2;
t = rightBC1(Amax);
fprintf('%e: function when Amax=%f\n',t,Amax)
A=Amin:.1:Amax;
for i=1:length(A)
xx(i) = A(i);
yy(i) = rightBC1(A(i));
end
figure(1)
hold on
plot(xx,yy),title('The function rightBC')
%pause
% now apply bisection to find the correct A value
tol = 1e-3;
Acorrect = bisect(@rightBC,Amin,Amax,tol);
fprintf('The correct slope at x=0 is %e\n',Acorrect)
z0 = [1; Acorrect];
[X,Z] = ode45(@ode,x,z0);
figure(2)
hold on
plot(X,Z(:,1),'r'),title('The solution y(x) with the correct A')
function r = rightBC(A)
% solve the system to get [y(x=2) when y'(0)=A]
% because the function is defined within the script file, the
% x vector is a global variable, and is hence available here.
% ode45 is a MATLAB routine to solve the system using
% the 4th order Runge-Kutta method
z0 = [1; A];
[X,Z] = ode45(@ode,x,z0);
% output: Z(:,1) is z1 (which is y) and Z(:,2) is z2 (which is y')
m = length(Z(:,1));
r = Z(m,1) - 6;
figure(2)
hold on
plot(X,Z(:,1)),title('The iterated solutions')
pause
end
function r = rightBC1(A)
% solve the system to get [y(x=2) when y'(0)=A]
% because the function is defined within the script file, the
% x vector is a global variable, and is hence available here.
% ode45 is a MATLAB routine to solve the system using
% the 4th order Runge-Kutta method
z0 = [1; A];
[X,Z] = ode45(@ode,x,z0);
% output: Z(:,1) is z1 (which is y) and Z(:,2) is z2 (which is y')
m = length(Z(:,1));
r = Z(m,1) - 6;
end
function dydx = ode(x,z)
dydx = [ z(2); x-4*z(2)+6*z(1)];
end
end