MATLAB lessons N.4¶
A simple implementation of Euler scheme (vector form): File Euler.m
:open:
1%
2% This routine approximate the solution of the IVP (Initial Value Problem)
3%
4% / y' = f(t,y)
5% |
6% | y(a) = y[0]
7% \
8%
9% PARAMETERS:
10% FUN = reference to user defined function
11% t0 = Starting integration time
12% y0 = Initial data of the IVP
13% h = time step
14% N = number of step to be performed
15%
16% OUTPUT
17% T = vector containing the nodes where
18% the approximate solution is computed
19% t = t0, t0+h, t0+3*h,..... t0+N*h
20% Y = matrix whose columsn containing the approximate solution at nodes T, i.e
21% Y(:,1) ~ y(t0),
22% Y(:,2) ~ y(t0+h),
23% Y(:,3) ~ y(t0+2*h),
24% .....
25% Y(:,N+1) ~ y(t0+N*h)
26%
27function [T,Y] = Euler( FUN, t0, y0, h, N )
28 % pre allcation of the output vector T and Y
29 T = zeros(N+1,1) ;
30 Y = zeros(size(y0,1),N+1) ; % size(y0,1) number of rows of y0
31
32 % set T(1) and Y(1) with the initial condition
33 T(1) = t0 ;
34 Y(:,1) = y0 ; % Y(:,1) access the first column of matrix Y
35
36 % loop performing the Euler Steps
37 for k=1:N
38 T(k+1) = T(k) + h ;
39 Y(:,k+1) = Y(:,k) + h * feval( FUN, T(k), Y(:,k) ) ;
40 % feval( FUN, T(k), Y(k) ) == FUN( T(k), Y(k) )
41 end
42
43end
A simple ODE to test the code: File problem1.m
:open:
1%
2% function for the IVP
3%
4% x'(t) = -y(t)
5% y'(t) = x(t)
6%
7% x(0) = 1
8% y(0) = 0
9%
10% / x \
11% Z = | |
12% \ y /
13%
14function f = problem1( t, Z )
15 f = zeros(2,1) ; % allocate for a 2 elements column vector
16 f(1) = -Z(2) ;
17 f(2) = Z(1) ;
18end
A simple implementation of Heun scheme (vector form): File Heun.m
:open:
1%
2% This routine approximate the solution of the IVP (Initial Value Problem)
3%
4% / y' = f(t,y)
5% |
6% | y(a) = y[0]
7% \
8%
9% PARAMETERS:
10% FUN = reference to user defined function
11% t0 = Starting integration time
12% y0 = Initial data of the IVP
13% h = time step
14% N = number of step to be performed
15%
16% OUTPUT
17% T = vector containing the nodes where
18% the approximate solution is computed
19% t = t0, t0+h, t0+3*h,..... t0+N*h
20% Y = matrix whose columsn containing the approximate solution at nodes T, i.e
21% Y(:,1) ~ y(t0),
22% Y(:,2) ~ y(t0+h),
23% Y(:,3) ~ y(t0+2*h),
24% .....
25% Y(:,N+1) ~ y(t0+N*h)
26%
27function [T,Y] = Heun( FUN, t0, y0, h, N )
28 % pre allcation of the output vector T and Y
29 T = zeros(N+1,1) ;
30 Y = zeros(size(y0,1),N+1) ; % size(y0,1) number of rows of y0
31
32 % set T(1) and Y(1) with the initial condition
33 T(1) = t0 ;
34 Y(:,1) = y0 ; % Y(:,1) access the first column of matrix Y
35
36 % loop performing the Euler Steps
37 for k=1:N
38 T(k+1) = T(k) + h ;
39 K1 = feval( FUN, T(k), Y(:,k) ) ;
40 ETA = Y(:,k) + h * K1 ; % Euler step y(k+1) = y(k) + h * f(t(k),y(k))
41 K2 = feval( FUN, T(k+1), ETA ) ;
42 Y(:,k+1) = Y(:,k) + (h/2) * (K1+K2) ; % Heun correction step
43 % y(k+1) = y(k) + (h/2) * (f(t(k),y(k))+f(t(k+1),y(k+1)))
44 end
45
46end
A simple implementation of Runge Kutta oder 4 (vector form): File RK4.m
:open:
1%
2% This routine approximate the solution of the IVP (Initial Value Problem)
3%
4% / y' = f(t,y)
5% |
6% | y(a) = y[0]
7% \
8%
9% PARAMETERS:
10% FUN = reference to user defined function
11% t0 = Starting integration time
12% y0 = Initial data of the IVP
13% h = time step
14% N = number of step to be performed
15%
16% OUTPUT
17% T = vector containing the nodes where
18% the approximate solution is computed
19% t = t0, t0+h, t0+3*h,..... t0+N*h
20% Y = matrix whose columsn containing the approximate solution at nodes T, i.e
21% Y(:,1) ~ y(t0),
22% Y(:,2) ~ y(t0+h),
23% Y(:,3) ~ y(t0+2*h),
24% .....
25% Y(:,N+1) ~ y(t0+N*h)
26%
27function [T,Y] = RK4( FUN, t0, y0, h, N )
28 % pre allcation of the output vector T and Y
29 T = zeros(N+1,1) ;
30 Y = zeros(size(y0,1),N+1) ; % size(y0,1) number of rows of y0
31
32 % set T(1) and Y(1) with the initial condition
33 T(1) = t0 ;
34 Y(:,1) = y0 ; % Y(:,1) access the first column of matrix Y
35
36 % loop performing the Euler Steps
37 for k=1:N
38 T(k+1) = T(k) + h ;
39 K1 = feval( FUN, T(k), Y(:,k) ) ;
40 K2 = feval( FUN, T(k)+h/2, Y(:,k)+(h/2)*K1 ) ;
41 K3 = feval( FUN, T(k)+h/2, Y(:,k)+(h/2)*K2 ) ;
42 K4 = feval( FUN, T(k)+h, Y(:,k)+h*K3 ) ;
43 Y(:,k+1) = Y(:,k) + (h/6) * (K1+2*(K2+K3)+K4) ;
44 end
45end
Exact solution to check the order: File exact1.m
:open:
1%
2% exact solution for the IVP
3%
4% x'(t) = -y(t)
5% y'(t) = x(t)
6%
7% x(0) = 1
8% y(0) = 0
9%
10% / x \
11% Z = | |
12% \ y /
13%
14% x(t) = cos(t)
15% y(t) = sin(t)
16%
17function Z = exact1( t )
18 Z = zeros(2,size(t,1)) ; % allocate for a 2 elements column vector
19 Z(1,:) = cos(t) ;
20 Z(2,:) = sin(t) ;
21end
A simple script to check the order: File checkError.m
:open:
1%
2% Compute approximate solution with a numerical scheme "METHOD"
3% and compute numerically the order
4%
5function checkError( METHOD )
6 %
7 h = 0.1 ;
8 N = 100 ;
9
10 % compute approximate solution
11 t0 = 0 ;
12 y0 = [1;0] ;
13 [T,Y] = feval ( METHOD, @problem1, t0, y0, h, N ) ;
14 Z = exact1( T ) ;
15 E1 = max(max(abs(Y-Z))) ;
16
17 % h <- h/2 N <- 2*N
18 h = h/2 ;
19 N = N*2 ;
20
21 [T,Y] = feval ( METHOD, @problem1, t0, y0, h, N ) ;
22 Z = exact1( T ) ;
23 E2 = max(max(abs(Y-Z))) ;
24
25
26 % h <- h/2 N <- 2*N
27 h = h/2 ;
28 N = N*2 ;
29
30 [T,Y] = feval ( METHOD, @problem1, t0, y0, h, N ) ;
31 Z = exact1( T ) ;
32 E3 = max(max(abs(Y-Z))) ;
33
34 fprintf( 1, 'log(E1/E2)/log(2) = %f\n', log(E1/E2)/log(2) ) ;
35 fprintf( 1, 'log(E2/E3)/log(2) = %f\n', log(E2/E3)/log(2) ) ;
36
37end