MATLAB lessons N.6

A simple implementation of Euler scheme (vector form): File BVPshoot.m

:open:
 1%
 2% Solve the BVP of the form
 3%
 4%   y'' + p y' + q y = r
 5%   y(a) = ya, y(b) = yb
 6%
 7% with N subinterval and the shooting method
 8%
 9%   h    = (b-a)/N
10%   x(k) = a + h*k,   
11%   y(k) ~ y(a+h*k),
12%
13%  p : reference of a function of one argument p(x)
14%  q : reference of a function of one argument q(x)
15%  r : reference of a function of one argument r(x)
16%
17%  BCa : vector of two element [ a , ya ] with left BC
18%  BCb : vector of two element [ b , yb ] with ight BC
19%
20%  N : number of subdivision of interval [a,b]
21%
22function [x,y] = BVPshoot( p, q, r, BCa, BCb, N )
23
24  global pFunction qFunction rFunction ; % store problem in global variables
25
26  pFunction = p ;
27  qFunction = q ;
28  rFunction = r ;
29
30  % extract the boundary conditions
31  a  = BCa(1) ;
32  b  = BCb(1) ;
33  ya = BCa(2) ;
34  yb = BCb(2) ;
35  % fill x
36  h  = (b-a)/N ;
37  x  = a + h*[0:N]' ;
38  % first problem y(a) = ya, y'(a) = 0
39  [x1,z] = ode45(@localOde1,x,[ya;0]) ;
40  
41  % second problem y(a) = 0, y'(a) = 1 and with r=0
42  [x2,w] = ode45(@localOde2,x,[0;1]) ;
43  
44  % lambda for the linear combination solution
45  lambda = (yb - z(end,1))/w(end,1) ; 
46  
47  % fill the y
48  y = z(:,1)+lambda*w(:,1) ;
49  
50end
51
52%
53%  Transformation of the second oder ODE
54%  
55%   y''(x) + p(x) y'(x) + q(x) y(x) = r(x)
56%  
57%  To a first order ODE (system)
58%
59%   y'(x) = v(x)
60%   v'(x) = r(x) - p(x) v(x) - q(x) y(x) 
61%
62%                    / y(x) \
63%  The vector Z(x) = |      |
64%                    \ v(x) /
65%
66function dZ = localOde1( x, Z )
67  global pFunction qFunction rFunction ;
68
69  dZ = [ Z(2) ; ...
70         feval(rFunction,x) ...
71         - feval(pFunction,x) *Z(2) ...
72         - feval(qFunction,x) * Z(1) ] ;
73
74end
75
76%
77%  Transformation of the second oder ODE
78%  
79%   y''(x) + p(x) y'(x) + q(x) y(x) = 0
80%  
81%  To a first order ODE (system)
82%
83%   y'(x) = v(x)
84%   v'(x) = - p(x) v(x) - q(x) y(x) 
85%
86%                    / y(x) \
87%  The vector Z(x) = |      |
88%                    \ v(x) /
89%
90function dZ = localOde2( x, Z )
91  global pFunction qFunction rFunction ;
92
93  dZ = [ Z(2) ; ...
94         - feval(pFunction,x) *Z(2) ...
95         - feval(qFunction,x) * Z(1) ] ;
96
97end

A simple BVP defined by the functions p1, q1, r1 to test the code: File p1.m

:open:
1function res = p1( x )
2  res = -x ;
3end

File q1.m

:open:
1function res = q1(x)
2  %res = -1 ;
3  res = -ones(size(x)) ;
4end

File r1.m

:open:
1function res = r1(x)
2  %res = 1 ;
3  res = ones(size(x)) ;
4end

Exact solution to check the order: File exact1.m

:open:
1function res = exact1(x)
2  C1  = (3*exp(-0.5) - 1) / erf( 1/sqrt(2)) ;
3  res = -1 + ( 1 + C1 * erf( x ./ sqrt(2) ) ) .* exp( (x.^2) ./ 2 ) ;
4end

A simple script to check the order: File checkError.m

:open:
 1BCa = [0;0] ;
 2BCb = [1;2] ;
 3
 4N = [ 10, 20, 40, 80, 160, 320, 640];%, 10000, 50000, 100000 ] ;
 5
 6for i=1:length(N)
 7  [x,y] = BVPshoot( @p1, @q1, @r1, BCa, BCb, N(i) );
 8  err(i) = norm( y - exact1( x ), inf ) ;
 9  disp( sprintf( 'error = %g\n', err(i) ) ) ;
10end