Applied Numerical Methods DN2221 - HS 2012

Computer Lab 2 Numerical Solution of Initial Value Problems

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Contents
1 2 3 Accuracy of Runge-Kutta Method Stability Investigation of a Runge-Kutta Method Parameter Studies of Solutions of an ODE System 2 3 7 9

A MATLAB code

A.1 Accuracy of a Runge-Kutta Method . . . . . . . . . . . . . . .

A.2 Constant Stepsize Experiment . . . . . . . . . . . . . . . . . . 12 A.3 Adaptive Stepsize Using MATLAB functions . . . . . . . . . . 15 A.4 Particle Flow Past a Cylinder . . . . . . . . . . . . . . . . . . 17 A.5 Motion of a Particle
B Bibliography

. . . . . . . . . . . . . . . . . . . . . . . 19
21

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Accuracy of Runge-Kutta Method

In the rst part of this lab we determine the accuracy of the given RungeKutta Method (RKM) in the lab description [1, p.21] in a numerical experiment, using the Van der Pol's dierential equation (VdPDE) [1, p.219] as system under investigation. First, we implement the given RKM to solve the VdPDE for constant stepsizes N = 10, 20, 40, 80, 160, 320 on the t-interval [0,1]. We then dene the local truncation error at t = 1 as eN = yN y (1), where y (1) yN =320 . Fig. 1 below shows a double logarithmic plot of the |eN | as a function of the (decreasing) stepsize h. Making a linear t to the logarithms of the calculated data in MATLAB, we nd that the order of the accuracy is approximately O(h3 ). This consistent with our expectations, since the applied RKM uses three coecients k1 , k2 and k3 . The MATLAB code for this part of the lab is included in the appendix (A.1).

Figure 1: Accuracy of the given RKM for VdPDE

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Stability Investigation of a Runge-Kutta Method

In the second part of this lab, we investigate the stability of the given RKM for Robertson's Problem given in [1, p.219]. First we use the RKM given in the rst part of the lab to solve Robertson's Problem with constant stepsizes corresponding to N = 125, 250, 500, 1000, 2000. Tab. 1 below shows the results of the calculations, where nan stands for not any number. In analogy to the discussion in [1, p.49], we conclude from tab. 1, that the smallest number of steps in order to achieve numerical stability is N = 1000. Fig. 2 below shows the solution trajectory for Robertson's Problem for N = 2000. Table 1: Results Constant Step-Size Experiment N 125 250 500 1000 2000 t

x1

x2

x3

1.0000    1.0000    1.0000    5 5 1.0000 0.96646 3.0746 10 3.0746 10 1.0000 0.96646 3.0746 105 3.0746 105

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Figure 2: Solution trajectory of Robertson's Problem for N = 2000

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Next, we perform an adaptive step-size experiment on Robertson's problem, using MATLAB functions. First, we use the non-sti IVP solver od23 on the t-interval [0,1] at dierent relative tolerances RelT ol = 103 , 104 , 105 , 106 and record the number of steps taken. Next, we run the sti IVP solver
ode23s

on the t-interval [0,1000] for the same tolerances and record the num-

ber of steps again. The recorded step-sizes for both experiments are summarized in tab. 2 below. Here we see, that the number of steps used in the rst experiment is an order of magnitude higher than the number of steps used in the second experiment. Moreover, the number of steps is practically constant for dierent relative tolerances in the rst experiment. On the contrary, the number of steps increases as the relative tolerance becomes smaller in the second experiment. Table 2: Results Constant Step-Size Experiment RelTol ode23 ode23s

103
866 30

104
867 37

105
868 48

106
868 61

Fig. A.3 below shows two graphs of the step-size h as a function of t at

RelT ol = 106 . The upper plot is generated using ode23 in the rst experiment, wheres the lower plot is generated using ode23s in the second experiment. The two plots visualize the observation of constant step-sizes (rst experiment) and varying step-sizes (second experiment) made above. The MATLAB code for this part of the lab is included in the appendix (A.2 and A.3)

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Figure 3: Comparison of step-sizes used by ode23 and ode23s to solve Robertson's problem.

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Parameter Studies of Solutions of an ODE System

In the last part of the lab, we perform a parameter study of the two problems given in the lab description [1, pp.220-221]. In the rst problem we calculate the ow curves of four dierent particles around a cylinder placed in an incompressible liquid on the t-interval [0,10]. In the second problem, we calculate the trajectories of a particle that is thrown up for dierent elevation angles and air resistance coecients. The MATLAB code for this part of the lab is included in the appendix (A.4 and A.5). Fig. 4 below shows the calculated ow curves for the rst problem, and g. 5 shows the calculated trajectories for the motion of the particle.

Figure 4: Parametric plot of ow curves of four dierent particles on t-interval [0,10]

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

Figure 5: Trajectories of particle for k = 0.02 (left) and k = 0.065 (right) for three dierent values of elevation angle = 30, 45 and 60 degrees.

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

A
A.1

MATLAB code
Accuracy of a Runge-Kutta Method

% f i r s t r e w r i t e 2nd o r d e r DE i n t o system o f f i r s t o r d e r DEs % d/ dt y_2 = y_1 e p s i l o n ( y_1^2 1) * y_2 % d/ dt y_1 = y_2 iterations = [10 ,20 ,40 ,80 ,160 ,320]; stepsize = 1./ i t e r a t i o n s ; data = z e r o s ( l e n g t h ( i t e r a t i o n s ) , 3 ) ; f 1 = i n l i n e ( ' y2 ' , ' y1 ' , ' y2 ' ) ; f 2 = i n l i n e ( ' y1 (y1^2 1 ) * y2 ' , ' y1 ' , ' y2 ' ) ; f o r k=1: l e n g t h ( i t e r a t i o n s ) h = stepsize (k) ; n = iterations (k) ; y1 = [ 1 : n + 1 ] ; y2 = [ 1 : n + 1 ] ; y1 ( 1 ) = 1 ; y2 ( 1 ) = 0 ; %implement runge k u t t a f o r system o f two c o u p l e d f i r s t o r d e r DEs f o r i =1:( n ) k1 = f 1 ( y1 ( i ) , y2 ( i ) ) ; m1 = f 2 ( y1 ( i ) , y2 ( i ) ) ; k2 = f 1 ( y1 ( i )+h * k1 , y2 ( i )+h *m1) ; %a r r a y t o s t o r e y1 %a r r a y t o s t o r e y2

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell m2 = f 2 ( y1 ( i )+h * k1 , y2 ( i )+h *m1) ;

10

k3 = f 1 ( y1 ( i ) +0.25 * h * k1 +0.25 * h * k2 , y2 ( i ) +0.25 * h * m1+0.25 * h *m2) ; m3 = f 2 ( y1 ( i ) +0.25 * h * k1 +0.25 * h * k2 , y2 ( i ) +0.25 * h * m1+0.25 * h *m2) ; y1 ( i +1) = y1 ( i ) + h / 6 * ( k1 + k2 + 4 * k3 ) ; y2 ( i +1) = y2 ( i ) + h / 6 * (m1 + m2 + 4 *m3) ; end %w r i t e r e s u l t s from i n n e r f o r l o o p i n t o data matrix data ( k , 2 ) = y1 ( n+1) ; data ( k , 3 ) = y2 ( n+1) ; c l e a r y1 ; c l e a r y2 ; end %c a l c u l a t e t h e e r r o r y_max = data ( l e n g t h ( i t e r a t i o n s ) , 2 ) ; e r r o r = abs ( data ( : , 2 ) y_max) ; error = error (1:5) ; stepsize = stepsize (1:5) ; %u s e data from l i n e a r f i t t o p l o t l i n e a r f i t i n t o same plot l i n f i t =l o g ( 4.3328) +3.0019 * l o g ( s t e p s i z e ) ; %p l o t e r r o r vs s t e p s i z e loglog ( stepsize , error , stepsize , l i n f i t ) t i t l e ( ' Accuracy RKM f o r VdPDE' ) xlabel ( ' stepsize h ') y l a b e l ( ' e r r o r i n comparison t o h = 1 / 3 2 0 ' ) ;

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

11

% %t h i s i s j u s t f o r e s t i m a t i n g t h e t h e e r r o r % s t e p s i z e=l o g ( s t e p s i z e ) % e r r o r=l o g ( e r r o r ) % plot ( stepsize , error )

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

12

% f i r s t r e w r i t e 2nd o r d e r DE i n t o system o f f i r s t o r d e r DEs % d/ dt y_2 = y_1 e p s i l o n ( y_1^2 1) * y_2 % d/ dt y_1 = y_2 %format l o n g ; c1 = 0 . 0 4 ; c2 = 1 0 ^ 4 ; c3 = 3 * 1 0 ^ 7 ; iterations = [125 ,250 ,500 ,1000 ,2000]; stepsize = 1./ i t e r a t i o n s ; data = z e r o s ( l e n g t h ( i t e r a t i o n s ) , 4 ) ; f 1 1 = i n l i n e ( ' 0.04 * y1 + ( 1 0 ^ 4 ) * y2 * y3 ' , ' y1 ' , ' y2 ' , ' y3 ' ) ; f 2 2 = i n l i n e ( ' 0 . 0 4 * y1 ( 1 0 ^ 4 ) * y2 * y3 ( 3 * 1 0 ^ 7 ) * ( y2 ) ^ 2 ' , ' y1 ' , ' y2 ' , ' y3 ' ) ; f 3 3 = i n l i n e ( ' 3 * 1 0 ^ 7 * ( y2 ) ^ 2 ' , ' y1 ' , ' y2 ' , ' y3 ' ) ; f o r k=1: l e n g t h ( i t e r a t i o n s ) h = stepsize (k) ; n = iterations (k) ; t = [ 1 : n+1]; y1 = [ 1 : n + 1 ] ; y2 = [ 1 : n + 1 ] ; y3 = [ 1 : n + 1 ] ; t (1) = 0; %u s e t a s parameter %a r r a y t o s t o r e y1 %a r r a y t o s t o r e y2

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell y1 ( 1 ) = 1 ; y2 ( 1 ) = 0 ; y3 ( 1 ) = 0 ; %implement runge k u t t a f o r system o f two c o u p l e d f i r s t o r d e r DEs f o r i =1:( n ) k1 = f 1 1 ( y1 ( i ) , y2 ( i ) , y3 ( i ) ) ; m1 = f 2 2 ( y1 ( i ) , y2 ( i ) , y3 ( i ) ) ; l 1 = f 3 3 ( y1 ( i ) , y2 ( i ) , y3 ( i ) ) ; k2 = f 1 1 ( y1 ( i )+h * k1 , y2 ( i )+h *m1, y3 ( i )+h * l 1 ) ; m2 = f 2 2 ( y1 ( i )+h * k1 , y2 ( i )+h *m1, y3 ( i )+h * l 1 ) ; l 2 = f 3 3 ( y1 ( i )+h * k1 , y2 ( i )+h *m1, y3 ( i )+h * l 1 ) ;

13

k3 = f 1 1 ( y1 ( i ) +0.25 * h * k1 +0.25 * h * k2 , y2 ( i ) +0.25 * h *m1+0.25 * h *m2, y3 ( i ) +0.25 * h * l 1 +0.25 * h * l 2 ) ; m3 = f 2 2 ( y1 ( i ) +0.25 * h * k1 +0.25 * h * k2 , y2 ( i ) +0.25 * h *m1+0.25 * h *m2, y3 ( i ) +0.25 * h * l 1 +0.25 * h * l 2 ) ; l 3 = f 3 3 ( y1 ( i ) +0.25 * h * k1 +0.25 * h * k2 , y2 ( i ) +0.25 * h *m1+0.25 * h *m2, y3 ( i ) +0.25 * h * l 1 +0.25 * h * l 2 ) ; t ( i +1) = t ( i ) + h ; y1 ( i +1) = y1 ( i ) + h / 6 * ( k1 + k2 + 4 * k3 ) ; y2 ( i +1) = y2 ( i ) + h / 6 * (m1 + m2 + 4 *m3) ; y3 ( i +1) = y3 ( i ) + h / 6 * ( l 1 + l 2 + 4 * l 3 ) ; end %w r i t e r e s u l t s from i n n e r f o r l o o p i n t o data matrix data ( k , 1 ) = t ( n+1) ; data ( k , 2 ) = y1 ( n+1) ; data ( k , 3 ) = y2 ( n+1) ; data ( k , 4 ) = y2 ( n+1) ;

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell %s t o r e s o l u t i o n t r a j e c t o r y f o r k=5o r n = 2000 t r a j 1 = y1 ; t r a j 2 = y2 ; t r a j 3 = y3 ; end %d i s p l a y r e s u l t s d i s p ( data ) %e x p o r t r e s u l t s t o l a t e x u s i n g l a t e x .m l i b r a r y from MathWorks . com l a t e x ( data , ' nomath ' ) loglog ( t , traj1 , t , traj2 , t , traj3 ) t i t l e ( ' s o l u t i o n t r a j e c t o r y R o b e r t s o n s Problem f o r N = 2000 ') a x i s ([10^( 4) 10^0 10^( 6) 1 0 ^ 1 ] ) x l a b e l ( ' parameter l n ( t ) ' ) y l a b e l ( ' s o l u t i o n components l n ( y ) ' ) ; l e g e n d ( 'A( t ) ' , 'B( t ) ' , 'C( t ) ' ) ; %p r i n t ( ' dpdf ' , ' r300 ' , ' . . / f i g u r e s / part_2a . jpg ' )

14

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

15

%c r e a t e column v e c t o r t o p l u g i n t o o d e 2 3 s format l o n g ; t o l e r a n c e = [ 1 e 3, 1 e 4, 1 e 5, 1 e 6 ] ; range = [ 0 1 ; 0 1 0 0 0 ] ;

%s o l v e t h e DE system with ode23 and r a n g e [ 0 1 ] f o r i =1:4 o p t i o n s = o d e s e t ( ' RelTol ' , t o l e r a n c e ( i ) ) ; [ T,Y] = ode23 ( @ r i g i d , [ 0 1 ] , [ 1 0 0 ] , o p t i o n s ) ; n_ode23_steps ( i ) = l e n g t h (T) 1; end % make an a r r a y f o r s t e p s i z e t o be p l o t t e d f o r j = 1 : ( l e n g t h (T) 1) h ( j ) = T( j +1)T( j ) ; time ( j ) = T( j +1) ; end subplot (2 ,1 ,1) p l o t ( time , h ) t i t l e ( ' R o b e r t s o n s Problem with ode23 and r a n g e [ 0 1] ') x l a b e l ( ' time t ' ) ylabel ( ' stepsize h ') %s o l v e t h e DE system with o d e 2 3 s and r a n g e [ 0 1 0 0 0 ] f o r i =1:4 o p t i o n s = o d e s e t ( ' RelTol ' , t o l e r a n c e ( i ) ) ;

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell [ T,Y] = o d e 2 3 s ( @ r i g i d , [ 0 1 0 0 0 ] , [ 1 0 0 ] , o p t i o n s ) ; n_ode23s_steps ( i ) = l e n g t h (T) 1; end % make an a r r a y f o r s t e p s i z e t o be p l o t t e d f o r j = 1 : ( l e n g t h (T) 1) h2 ( j ) = T( j +1)T( j ) ; time2 ( j ) = T( j +1) ; end subplot (2 ,1 ,2) p l o t ( time2 , h2 )

16

t i t l e ( ' R o b e r t s o n s Problem with o d e 2 3 s and r a n g e [ 0 1000] ') x l a b e l ( ' time t ' ) ylabel ( ' stepsize h ') f u n c t i o n [ dy ] = r i g i d ( t , y ) k1 = 0 . 0 4 ; k2 = 1 0 ^ 4 ; k3 = 3 * 1 0 ^ 7 ; dy = z e r o s ( 3 , 1 ) ; dy ( 1 ) = k1 * y ( 1 ) + k2 * y ( 2 ) * y ( 3 ) ; dy ( 2 ) = k1 * y ( 1 ) k2 * y ( 2 ) * y ( 3 ) k3 * ( y ( 2 ) ) ^ 2 ; dy ( 3 ) = k3 * ( y ( 2 ) ) ^ 2 ; end

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

17

N=1000; h=10/N; y l i s t = [0.2 0.6 1 1 . 6 ] ; f o r k =1:4 u = [ 4; y l i s t ( k ) ] ; r e s u l t =[u ] ; h l i s t = [ 0 ] ; f o r j =1:N; u1 = [ u ( 1 ) ; u ( 2 ) ] ; k1 = C y l i n d e r ( u1 ) ; u2 = [ u ( 1 )+h * k1 ( 1 ) ; u ( 2 )+h * k1 ( 2 ) ] ; k2 = C y l i n d e r ( u2 ) ; u3 = [ u ( 1 ) +(h / 4 ) * k1 ( 1 ) +(h / 4 ) * k2 ( 1 ) ; u ( 2 ) +(h / 4 ) * k1 ( 2 ) +(h / 4 ) * k2 ( 2 ) ] ; k3 = C y l i n d e r ( u3 ) ; u = u+(h / 6 ) * ( k1+k2+4* k3 ) ; r e s u l t =[ r e s u l t u ] ; h l i s t =[ h l i s t j * h ] ; end h o l d on plot ( r e s u l t ( 1 , : ) , r e s u l t ( 2 , : ) ) ; axis equal ;

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell t i t l e ( ' p a r a m e t r i c p l o t o f p a r t i c l e f l o w around cylinder ' ) ; xlabel ( 'x ' ) ; ylabel ( 'y ' ) ; end

18

f u n c t i o n [ u_prime ] = C y l i n d e r ( u ) R = 2; r_sq = ( u ( 1 ) ^2+u ( 2 ) ^2) ^ 2 ; u_prime = [1 R^2 * ( u ( 1 ) ^2u ( 2 ) ^2) / r_sq ; 2*R^2 * u ( 1 ) * u ( 2 ) / r_sq ] ; end

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

19

h=10^( 3) ; f o r i =1:2; k =[0.02 0 . 0 6 5 ] ; f o r j =1:3 a l p h a = [ p i /6 p i /4 p i / 3 ] ; r = [0; 1.5]; v = [ 2 0 * c o s ( a l p h a ( j ) ) ; 20 * s i n ( a l p h a ( j ) ) ] ; r e s u l t =[ r ] ; h l i s t = [ 0 ] ; w h i l e r ( 2 )>0 v1 = [ v ( 1 ) ; v ( 2 ) ] ; k1 = Motion ( v1 , k ( i ) ) ; v2 = [ v ( 1 )+h * k1 ( 1 ) ; v ( 2 )+h * k1 ( 2 ) ] ; k2 = Motion ( v2 , k ( i ) ) ; v3 = [ v ( 1 ) +(h / 4 ) * k1 ( 1 ) +(h / 4 ) * k2 ( 1 ) ; v ( 2 ) +(h / 4 ) * k1 ( 2 ) +(h / 4 ) * k2 ( 2 ) ] ; k3 = Motion ( v3 , k ( i ) ) ; v = v+(h / 6 ) * ( k1+k2+4* k3 ) ; r = r + h*v ; result = [ result r ] ;

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell end h o l d on subplot (1 ,2 , i ) ; plot ( result (1 ,:) , result (2 ,:) ) ; t i t l e ( [ ' k = ' , num2str ( k ( i ) ) ] ) ; x l a b e l ( ' x ' ) ; ylabel ( 'y ') ; a x i s ( [ 0 30 0 1 5 ] ) ;

20

end end

f u n c t i o n [ v_prime ] = Motion ( v , k ) r_prime = ( v ( 1 ) ^2+v ( 2 ) ^2) ^ ( 1 / 2 ) ; v_prime = [ k * v ( 1 ) * r_prime ; 9.81 k * abs ( v ( 2 ) ) * r_prime ] ; end

DN2221 - HS 2012 - Computer Lab 2 - S.Arnoldt & J. Norell

21

Bibliography

References
[1] L. Edsberg, 2008. Introduction to Computation and Modelling for Differential Equations.

1st edition. New Jersey: John Wiley & Sons.

[2] MATLAB, 2012. MATLAB R2012a (7.14.0.739). Natick, MA: The MathWorks Inc.