The Shooting Method for Boundaryvalue Problems Douglas Wilhelm

The Shooting Method for Boundary-value Problems Douglas Wilhelm Harder, M. Math. LEL Department of Electrical and Computer Engineering University of Waterloo, Ontario, Canada ece. uwaterloo. ca dwharder@alumni. uwaterloo. ca © 2012 by Douglas Wilhelm Harder. Some rights reserved.

The Shooting Method for Boundary-value Problems Outline This topic to approximating boundary-value problems – We will describe boundary-value problems (BVPs) – We will look at solutions with linear ordinary differential equations (ODEs) – We will consider solutions for non-linear ODEs • This will require successive approximations using the secant method 2

The Shooting Method for Boundary-value Problems Outcomes Based Learning Objectives By the end of this laboratory, you will understand: – Boundary-value problems (BVPs) – How to use an initial-value problem (IVP) solver to approximate solutions to BVPs – The solution to linear BVPs – The application of the secant method to solve non-linear BVPs 3

The Shooting Method for Boundary-value Problems Using Matlab These methods assume that the student has written the initial-value problem solver pd 45 with the signature dp 45( f, x_rng, u 0, h, eps_abs ) which uses the Dormand-Prince method If this routine is not available, you are welcome to use the built-in Matlab routine ode 45 which has the similar signature ode 45( f, x_rng, u 0 ) 4

The Shooting Method for Boundary-value Problems 2 nd-order ODEs A boundary-value problem in one dimension is any 2 ndorder ODE F(x, u(x), u(1)(x), u(2)(x)) = 0 with two constraints u(a) = ua u(b) = ub In general, we will look at functions of the form u(2)(x) = f(x, u(x), u(1)(x)) 5

The Shooting Method for Boundary-value Problems 2 nd-order ODEs Consider the corresponding initial-value problem for this 2 nd-order ODE F(x, u(x), u(1)(x), u(2)(x)) = 0 with two constraints u(a) = ua u(1)(a) = ua(1) Here we specify the slope at the left-hand point x = a 6

The Shooting Method for Boundary-value Problems 2 nd-order ODEs Thus, a boundary-value problem could be restated as: – Given one initial condition u(a) = ua, what slope is required at that initial point so that the solution to that initial condition passes through the point (b, ub)? (b, ub) (a, ua) 7

The Shooting Method for Boundary-value Problems 2 nd-order ODEs We will consider two possible cases: – When the ODE is linear, and – When it is not (b, ub) (a, ua) 8

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs A 2 nd-order ODE is said to be linear if it can be written as u(2)(x) + q(x) u(1)(x) + r(x) u(x) = g(x) or u(2)(x) = g(x) – q(x) u(1)(x) – r(x) u(x) We will also consider the corresponding homogenous LODE: u(2)(x) + q(x) u(1)(x) + r(x) u(x) = 0 or u(2)(x) = –q(x) u(1)(x) – r(x) u(x) 9

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs To demonstrate 2 nd-order LODEs, we will consider the example The corresponding homogeneous LODE is the 10

The Shooting Method for Boundary-value Problems 11 2 nd-order Linear ODEs Recall that if ug(x) is a solution to a 2 nd-order LODE and u 0(x) is a solution to the corresponding homogenous LODE, then ug(x) + c·u 0(x) is also a solution the LODE g(x) 0

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs Suppose we have redefine a boundary-value problem in terms of two IVPs, the first on the LODE with constraints and the second on the homogenous LODE with constraints Let ug(x) and u 0(x) be the solutions, respectively 12

The Shooting Method for Boundary-value Problems 13 2 nd-order Linear ODEs For example, given the LODE from before, we create the two IVPs (b, ub) (a, ua) The plots of the solutions: ug(x) u 0(x)

The Shooting Method for Boundary-value Problems 14 2 nd-order Linear ODEs Having found these two solutions, ug(x) and u 0(x), we now want to find a linear combination of these two such that it satisfies the second boundary value (b, u ) b – For example, here we see ug(x), ug(x) + u 0(x), and ug(x) – u 0(x) It seems ug(x) + 2 u 0(x) will pass close to the second boundary point ug(x) + u 0(x) (a, ua) ug(x) – u 0(x)

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs We want to find the appropriate linear combination to match our second boundary condition First, we note We want Solving this for c yields 15

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs Easy enough; however, recall that we must find ug(x) and u 0(x) but if we use dp 45, the step sizes may be different… Solution: define a system of two uncoupled equations and unknowns and solve them simultaneously u(2)(x) = g(x) – q(x) u(1)(x) – r(x) u(2)(x) = 0 – q(x) u(1)(x) – r(x) u(x) 16

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs Using our techniques from Laboratory 5, we can write this as a system of IVPs: w 1(1)(x) = w 2 (x) w 2(1)(x) = g(x) – q(x) w 2(x) – r(x) w 1(x) w 3(1)(x) = w 4 (x) w 4(1)(x) = 0 – q(x) w 4(x) – r(x) w 3(x) w 1(a) = ua w 2(a) = 0 w 3(a) = 0 w 4(a) = 1 ug(x) = w 1(x) ug(1)(x) = w 2(x) u 0(x) = w 3(x) u 0(1)(x) = w 4(x) 17

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs For example, consider the boundary-value problem The Matlab function would be: function [dw] = f 6 a( x, w ) dw = [w(2); 1 - sin(x)*w(2) - w(1); w(4); 0 - sin(x)*w(4) - w(3)]; end with the initial conditions [u_a, 0, 0, 1]' 18

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs For our example, we would run the following code: a = 2; b = 3; u_a = 1. 5; u_b = 2. 5; [x 6 a, u 6 a] = dp 45( @f 6 a, [a, b], [u_a 0 0 1]', 0. 01, 1 e-8 ); c 6 a = (u_b - u 6 a(1, end))/u 6 a(3, end) c 6 a = 1. 991628911663175 u 6 a_soln = u 6 a(1, : ) + c 6 a*u 6 a(3, : ); plot( x 6 a, u 6 a_soln, 'r' ) ylim( [0, 2. 5] ) u 6 a_soln(1) ans = 1. 50000000 u 6 a_soln(end) ans = 2. 50000000 19

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs Maple can find an answer, but it’s not nice… > a : = 2; b : = 3; u_a : = 1. 5; u_b : = 2. 5; > dsolve( {(D@@2)(u)(x) + sin(x)*D(u)(x) + u(x) = 1, u(a)=u_a, u(b)=u_b} ); 20

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs It takes Maple about two minute to plot the same figure – Note: this is not an issue with Maple—it is giving you the exact solution and not a numeric approximation 21

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs: Example For example, consider the boundary-value problem u(2)(x) + 4 u(1)(x) + 7 x u(x) = sin(x) u(2) = 1. 5 u(3) = 2. 5 The Matlab function would be: function [dw] = f 6 b( x, w ) dw = [w(2); sin(x) - 4*w(2) - 7*x*w(1); w(4); 0 - 4*w(4) - 7*x*w(3)]; end 22

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs: Example For example, consider a = 2; b = [x 6 b, u 6 b] u 6 b_soln = plot( x 6 b, grid 3; u_a = 1. 5; u_b = 2. 5; = dp 45( @f 6 b, [a, b], [u_a 0 0 1]', 0. 01, 1 e-4 ); u 6 b(1, : ) + (u_b - u 6 b(1, end))/u 6 b(3, end)*u 6 b(3, : ); u 6 b_soln, 'o-' ); 23

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs: Example Maple can find an answer, but it’s not nice… > a : = 2; b : = 3; u_a : = 1. 5; u_b : = 2. 5; > dsolve( {(D@@2)(u)(x) + 4*D(u)(x) + 7*x*u(x) = sin(x), u(a) = u_a, u(b) = u_b}, u(x) ); 24

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs: Example Comparing the results with Maple’s plot, we see they are similar 25

The Shooting Method for Boundary-value Problems 2 nd-order Linear ODEs: Example If we repeat using b = 4, we get significantly more varying, but still similar, graphs 26

The Shooting Method for Boundary-value Problems 2 nd-order Non-linear ODEs What happens if the ODE is not linear? – In this case, there is no guarantee that we can find a solution using this simple technique – Instead, we will use u(2)(x) = f(x, u(x), u(1)(x)) and let us(x) be the solution to the IVP u(a) = ua u(1) (a) = s Now, define a new function err(s) = us(b) – ub 27

The Shooting Method for Boundary-value Problems 2 nd-order Non-linear ODEs When this function has a zero, we have a solution to the boundary-value problem err(s) = us(b) – ub What is err(s) in Matlab? function [du_b] = err_shot( s ) [t_out, u_out] =. . . dp 45( @f, [a, b], [u_a, s], h, eps_abs ); du_b = u_out(1, end) - u_b; end 28

The Shooting Method for Boundary-value Problems 29 2 nd-order Non-linear ODEs When this function has a zero, we have a solution to the boundary-value problem err(s) = us(b) – ub In reality, however, it will require a number of other parameters: function [du_b] = err_shot(s, f, x_rng, u_bndry, h, eps_abs) u_a = u_bndry(1); u_b = u_bndry(2); [~, u_out] = dp 45( f, x_rng, [u_a, s]', h, eps_abs ); du_b = u_out(1, end) - u_b; end

The Shooting Method for Boundary-value Problems Secant Method We will use the secant method for fining this root: – Suppose we are trying to find a root of a real-valued function of a real variable err(s) – Suppose s 1 and s 2 are two initial approximations of the root err(s) s 2 s 1 30

The Shooting Method for Boundary-value Problems Secant Method We can find a better approximation of the root by finding the interpolating straight line that passes through these points err(s) s 2 s 1 31

The Shooting Method for Boundary-value Problems Secant Method The formula for this new point is err(s) s 3 s 2 s 1 32

The Shooting Method for Boundary-value Problems Secant Method In the first case, we considered linear ODEs – Essentially, what we are doing here is asking: Given these two approximations, what would our next approximation be if the ODE was linear? err(s) s 3 s 2 s 1 33

The Shooting Method for Boundary-value Problems Secant Method Now, suppose we have s 2 and s 3, the next step is to use these approximations to find the next approximation, s 4, and so on. . . err(s) s 3 s 2 s 1 34

The Shooting Method for Boundary-value Problems Secant Method Thus, we will continue to iterate, using the secant method at each step. . . With any iterative numerical method, we must know under what conditions we will: – Halt with success, and – Indicate a failure to find a solution After all, not all iterative numerical methods converge. . . 35

The Shooting Method for Boundary-value Problems Secant Method For the secant method, we must have two approximations, s 1 and s 2 First, if |err(s 1)| < |err(s 2)|, we will swap s 1 and s 2 – We’re assuming that s 2 is more accurate 36

The Shooting Method for Boundary-value Problems Secant Method We will iterate at most Nmax times: – With each iteration, we will approximate a new point s – If |s 2 – s| < estep and |err(s)| < eabs, we will assume we are finished and we will return the value s (shooting, however, will take one more step) – Otherwise, we will set s 1 = s 2 and s 2 = s If we have iterated Nmax times and not found a solution, we will throw an exception with an appropriate message 37

The Shooting Method for Boundary-value Problems The Shooting Method Normally, the secant method simply returns the slope Your function, shooting, however, will have to return the x- and u-values corresponding to that optimal slope – Thus, your last steps before you successfully return will be: [x_out, u_out] = dp 45( f, x_rng, [u_bndry(1), s]', h, eps_abs ); return; 38

The Shooting Method for Boundary-value Problems The Shooting Method The arguments of function [x_out, u_out] = shooting( s 1, s 2, f, x_rng, u_bndry, . . . h, eps_abs, eps_step, N_max ) are as follows: s 1 and s 2 are the two approximations of the slopes f is a function handle for the differential equation: u(2)(x) = f(x, u(x), u(1)(x)) [a, b] is a row vector defining the range on which we are approximating the boundary-value problem [ua, ub] is a row vector defining the boundary conditions: u(a) = ua and u(b) = ub h is the initial step size for the function dp 45 eabs is the parameter passed to dp 45 and is also the parameter used by the secant method estep is the parameter used by the secant method for the step size Nmax is the parameter used by the secant method to define the maximum number of iterations 39

The Shooting Method for Boundary-value Problems The Shooting Method Question: what do we use as initial conditions? Answer: as an engineer, you will understand your problem; you should know the expected behaviour around the point x = a – In these laboratories, you will be given the initial points – In reality, you will understand what they should be when you get there 40

The Shooting Method for Boundary-value Problems The Shooting Method Consider the following non-linear boundary-value problem: The corresponding function is function [dw] = f 6 c( x, w ) dw = [w(2); sin(x) - 4*w(2)*w(1) + 2*x*w(1)^2]; end 41

The Shooting Method for Boundary-value Problems The Shooting Method When I run the function [x 6 c, u 6 c] = shooting( -3, -3. 1, @f 6 c, [2, 3], [1. 5, 2. 5]', . . . 0. 01, 1 e-6, 20 ); I get a solution after four iterations where the approximations of the initial slopes are: – 3. 326745004558498 – 3. 323136496135223 – 3. 323176635129875 – 3. 323176642445732 These last two slopes are sufficiently close enough and the u 6 c(1, end) = 2. 499999991 close enough to 2. 5 42

The Shooting Method for Boundary-value Problems The Shooting Method We can look at the plot and the result: plot( xs, us(1, : ) ) us(1, end) ans = 2. 499999991 length( xs ) ans = 36 43

The Shooting Method for Boundary-value Problems The Shooting Method You might get slightly different values depending on your implementation of the secant method; however, the final result should be reasonably close 44

The Shooting Method for Boundary-value Problems 45 The Shooting Method Plotting approximations with the various slopes: slopes = [-3 -3. 1 -3. 326745004558498 -3. 323136496135223. . . -3. 323176635129875 -3. 323176642445732 ]; for s = slopes [x 6 s, u 6 s] = dp 45( @f 6 c, [2, 3], [1. 5, s]', 0. 01, 1 e-6 ); plot( x 6 s, u 6 s(1, : ) ) u 6 s(1, end) end s 1 = 3 s 2 = 3. 1 All others. . .

The Shooting Method for Boundary-value Problems The Shooting Method If you plot the relative errors of the approximation of u 6 s(1, end) as an approximation of ub = 2. 5, we observe the fast (i. e. , O(hf)) convergence of the secant method s u 6 s(1, end) Absolute Error – 3 2. 668638013935938 1. 686 × 10– 1 – 3. 1 2. 617026509067240 1. 170 × 10– 1 – 3. 326745004558498 2. 498107476136670 1. 893 × 10– 3 – 3. 323136496135223 2. 500021288156047 2. 129 × 10– 5 – 3. 323176635129875 2. 500000003879338 3. 879 × 10– 9 – 3. 323176642445732 2. 499999991 7. 105 × 10– 15 46

The Shooting Method for Boundary-value Problems Hints You will work with various values of s and err(s) – Instead of continually recalculating value err(s), just calculate err 1 = err_shot( s 1 ); err 2 = err_shot( s 2 ); and when you calculate a new value of s, just use a simpler expression with the variables s 1, s 2, err 1 and err 2 and then immediately calculate errs = err( s ); – When you update s 1 and s 2, update err 1 and err 2 47

The Shooting Method for Boundary-value Problems Summary We have looked at using the function dp 45 to approximate a boundary-value problem – The shooting method converts a BVP into an IVP – For linear IVPs, we use the ODE and the homogeneous ODE and find a linear combination of the two solutions – For non-linear IVPs, we use dp 45 and the secant method to guide us to an approximation 48

The Shooting Method for Boundary-value Problems References [1] Glyn James, Modern Engineering Mathematics, 4 th Ed. , Prentice Hall, 2007. [2] Glyn James, Advanced Modern Engineering Mathematics, 4 th Ed. , Prentice Hall, 2011. [3] John H. Mathews and Kurtis D. Fink, Numerical Methods using Matlab, 4 th Ed. , Prentice Hall, 2004, pp. 529 -534. [4] Shooting Method, http: //en. wikipedia. org/wiki/Shooting_method. 49