Solving a 2nd order ODE with ode45

Contents

Initial value problem

We consider an initial value problem for a 2nd order ODE:

$$y''-y'+3y=t,\qquad y(0)=1,\quad y'(0)=-2$$

and we want to find the solution y(t) for t in [0,4].

We first have to rewrite this as a 1st order system: Let $y_1(t):=y(t)$ and $y_2(t)=y'(t)$, then we obtain

$$\left[\begin{array}{c}y_1' \\ y_2'\end{array}\right] = \left[\begin{array}{c}y_2 \\ t+y_2-3y_1\end{array}\right], \qquad \left[\begin{array}{c}y_1(0) \\ y_2(0)\end{array}\right] = \left[\begin{array}{c} 1\\ -2\end{array}\right]$$

Now we can define a vector valued function f(t,y) and an initial vector y0.

f = @(t,y) [y(2); t+y(2)-3*y(1)]; % define function f(t,y)
t0 = 0;  y0 = [1;-2];             % initial condition with vector y0

Solve IVP with ode45

ode45 returns a vector ts of t values and an array ys: each row of ys contains the values for $y_1$ and $y_2$.

Note that ys(end,1) is the approximation for $y_1(T)=y(T)$ and ys(end,2) is the approximation for $y_2(T)=y'(T)$.

ys(:,1) gives the values for $y_1(t)=y(t)$, ys(:,1) gives the values for $y_2(t)=y'(t)$. Therefore we can plot the function $y(t)$ with plot(ts,ys(:,1)

T = 4;                            % final t-value
[ts,ys] = ode45(f,[t0,T],y0);     % solve IVP

fprintf('y(T) = %g, y''(T) = %g\n',ys(end,1),ys(end,2))

disp('       t           y1          y2')
disp([ts,ys])                     % table of t and y values of solution

plot(ts,ys(:,1),'b')              % plot solution y(t)
title('Solution y(t) of IVP')
xlabel('t'); grid on

y(T) = 3.36902, y'(T) = -21.7257
       t           y1          y2
            0            1           -2
     0.020095       0.9588         -2.1
      0.04019      0.91561      -2.1992
     0.060285      0.87043      -2.2972
      0.08038      0.82329       -2.394
      0.18038      0.56069      -2.8524
      0.28038      0.25469      -3.2587
      0.38038    -0.088704      -3.5947
      0.48038     -0.46145      -3.8424
      0.58038     -0.85374       -3.985
      0.68038      -1.2543      -4.0065
      0.78038      -1.6505      -3.8935
      0.88038      -2.0283      -3.6356
      0.98038      -2.3727       -3.226
       1.0804      -2.6683      -2.6612
       1.1804      -2.8997      -1.9424
       1.2804      -3.0518      -1.0763
       1.3804      -3.1107    -0.075287
       1.4804      -3.0632        1.044
       1.5804      -2.8986       2.2589
       1.6804      -2.6089       3.5411
       1.7804      -2.1894       4.8569
       1.8804       -1.638       6.1687
       1.9804      -0.9571       7.4349
       2.0804     -0.15369       8.6109
       2.1804      0.76053       9.6502
       2.2804       1.7699       10.505
       2.3804       2.8539        11.13
       2.4804       3.9871       11.481
       2.5804       5.1398       11.518
       2.6804       6.2789       11.207
       2.7804       7.3685        10.52
       2.8804         8.37        9.441
       2.9804       9.2438       7.9625
       3.0804       9.9496        6.088
       3.1804       10.449       3.8344
       3.2804       10.705       1.2335
       3.3804       10.686      -1.6683
       3.4804       10.364      -4.8123
       3.5804       9.7173      -8.1258
       3.6804       8.7349      -11.522
       3.7603       7.7061      -14.229
       3.8402       6.4631       -16.87
       3.9201       5.0132       -19.39
            4        3.369      -21.726


Published with MATLAB® R2017a