% start off with an example for solving an ODE BVP by shooting. % WARNING: this example has a basis flaw; see below. % BVP via shooting % y'' = -y on [a .. b] = [0 .. pi], y(0) = 0, y(pi) = 1 % convert to first-order system: % z := [y;y'] >> f = inline('[z(2); - z(1)]','t','z') f = Inline function: f(t,z) = [z(2); - z(1)] >> ya = 0; yb = 1; a = 0; b = pi; % treat the value v at b of the (first component of the) solution as a % function v = v(s) of the slope s (i.e., the second component) at a % Choose a starting guess for the initial slope >> s1 = -1; >> [tsol,zsol1] = ode45(f,[a b],[ya;s1]); v1 = zsol1(end,1) v1 = 1.2628e-006 % a second guess: >> s2 = -100; >> [tsol,zsol2] = ode45(f,[a b],[ya;s2]); v2 = zsol2(end,1) v2 = 1.2628e-004 % we have two values ; ready to use the Secant method: >> s3 = s2 - (v2-1)/(v2 - v1)*(s2-s1) s3 = -7.9192e+005 % use this new value: >> [tsol,zsol3] = ode45(f,[a b],[ya;s3]); v3 = zsol3(end,1) v3 = 1.0062 %% well, we are almost there: >> s4 = s3 - (v3-1)/(v3 - v2)*(s3-s2) s4 = -7.8702e+005 >> [tsol,zsol4] = ode45(f,[a b],[ya;s4]); v4 = zsol4(end,1) v4 = 1.0032 % We might have hoped for a better results here, % ... but now consider the actual numbers involved here: a slope of -8e05 % seems remarkably peculiar! % Look at the mathematical problem we start off with. We are looking for a % solution to the homogeneous problem y'' + y = 0 , wishing to % pick a particular solution by imposing two side conditions on the general % solution. You have learned in the prereqs for this course that the general % solution to this particular problem is % c1*sin(t) + c2*cos(t) % Imposing the condition y(0) = 0 implies that c2 = 0 , leaving us with % c1*sin(t) % and this function vanishes also at t = pi. In other words, our original % problem HAS NO SOLUTION. No wonder the numerical problem is struggling so % hard. % We have here a specific example of the general fact that even a linear % ODE has solutions only locally, i.e., on a small enough interval. In this % particular case, if we switch to the smaller interval [a .. b] = [0 .. pi/2], % then the two-point boundary value problem has a unique solution, as we can % check by looking at the corresponding linear system % % c1*sin(a) + c2*cos(a) = ya % c1*sin(b) + c2*cos(b) = yb % % For a = 0, b = pi, the coefficient matrix is [0 1; 0 -1], and we can % see that it is not invertible. % For a = 0, b = pi/2 , the coefficient matrix is [0 1; 1 0] , and that % matrix is trivially invertible. The solution is c1 = yb, c2 = ya .