Numerical solution of IVP for 2nd order ODE
We want to find the solution of the IVP
for
. This is a second order ODE.
We first have to rewrite this as a first order system:
Let
with the vector valued function
. Note that We obtain the first order system
where we have
We define the function
in Matlab: f = @(t,y) [ y(2) ; t-3*y(2)-2*y(1) ];
Our initial vector is
Perform one step of Euler method with 
We evaluate
: and obtain for
the approximation 
We obtain the approximations
. Use 10 steps of the Euler method for 
Here the step size is
. The function EulerMethod is defined at the end of this file.
[ts,ys] = EulerMethod(f,[0,5],y0,10);
disp([ts,ys])
0 2.0000 -1.0000
0.5000 1.5000 -1.5000
1.0000 0.7500 -0.5000
1.5000 0.5000 0
2.0000 0.5000 0.2500
2.5000 0.6250 0.3750
3.0000 0.8125 0.4375
3.5000 1.0312 0.4688
4.0000 1.2656 0.4844
4.5000 1.5078 0.4922
5.0000 1.7539 0.4961
Note that ys has two columns containing the values of
and
. Plot the Euler approximation for
: xlabel('t'); title('Euler approximation for IVP');hold on
Using ode45
Find the solution
for
: [ts,ys] = ode45(f,[0,5],y0);
disp([ts,ys]) % shows t, y(t), y'(t) in each row (scroll to see all rows)
0 2.0000 -1.0000
0.0502 1.9486 -1.0430
0.1005 1.8954 -1.0728
0.1507 1.8410 -1.0910
0.2010 1.7860 -1.0992
0.2876 1.6909 -1.0939
0.3742 1.5972 -1.0688
0.4608 1.5062 -1.0285
0.5474 1.4192 -0.9772
0.6394 1.3322 -0.9145
0.7314 1.2512 -0.8461
0.8234 1.1766 -0.7741
0.9154 1.1087 -0.7006
1.0188 1.0406 -0.6183
1.1222 0.9809 -0.5374
1.2256 0.9294 -0.4589
1.3289 0.8858 -0.3837
1.4410 0.8472 -0.3067
1.5531 0.8169 -0.2344
1.6653 0.7945 -0.1671
1.7774 0.7792 -0.1048
1.8955 0.7705 -0.0445
2.0136 0.7685 0.0105
2.1317 0.7728 0.0607
2.2499 0.7827 0.1062
2.3612 0.7967 0.1450
2.4726 0.8149 0.1804
2.5840 0.8368 0.2124
2.6954 0.8620 0.2413
2.8204 0.8941 0.2706
2.9454 0.9296 0.2966
3.0704 0.9681 0.3198
3.1954 1.0094 0.3404
3.3204 1.0531 0.3587
3.4454 1.0990 0.3750
3.5704 1.1468 0.3894
3.6954 1.1963 0.4022
3.8204 1.2473 0.4135
3.9454 1.2996 0.4236
4.0704 1.3531 0.4324
4.1954 1.4077 0.4403
4.3204 1.4631 0.4473
4.4454 1.5194 0.4534
4.5704 1.5765 0.4589
4.6954 1.6341 0.4637
4.7715 1.6695 0.4663
4.8477 1.7052 0.4688
4.9238 1.7409 0.4710
5.0000 1.7769 0.4732
Note that ys has two columns containing the values of
and
. Show the value for
: Plot the solution
: plot(ts,ys(:,1),'k'); hold off
title('approximations from ode45 (black) and Euler method (blue)')
Code for Euler method
We modify the code from the scalar case so that it works for systems: for each time we store the vector
as a row of the array ys. function [ts,ys] = EulerMethod(f,interval,y,n)
% [ts,ys] = EulerMethod(f,[t0,T],y0,n)
% initial value problem y'=f(t,y), y(t0) = y0
% use n steps of Euler method to approximate solution in interval [t0,T]
t = interval(1); T = interval(2);
m = size(y,1); % size of vector y
ts = zeros(n+1,1); ys = zeros(n+1,m); % initialize vector for t-values, array for y-values
ts(1) = t; ys(1,:) = y.'; % y.' is transpose, gives row vector
ts(k+1) = t; ys(k+1,:) = y.'; % store t,y in vector ts, array ys