MA342 Handout: Solve system of ODEs Using Octave function LSODE =========================================================================== [X, ISTATE, MSG] = lsode (FCN, X_0, T, T_CRIT) Solve the set of differential equations dx -- = f(x, t) dt with x(t_0) = x_0 The solution is returned in the matrix X, with each row corresponding to an element of the vector T. The first element of T should be t_0 and should correspond to the initial state of the system X_0, so that the first row of the output is X_0. The first argument, FCN, is a string, inline, or function handle that names the function f to call to compute the vector of right hand sides for the set of equations. The function must have the form XDOT = f (X, T) in which XDOT and X are vectors and T is a scalar. If FCN is a two-element string array or a two-element cell array of strings, inline functions, or function handles, the first element names the function f described above, and the second element names a function to compute the Jacobian of f. The Jacobian function must have the form JAC = j (X, T) in which JAC is the matrix of partial derivatives | df_1 df_1 df_1 | | ---- ---- ... ---- | | dx_1 dx_2 dx_N | | | | df_2 df_2 df_2 | | ---- ---- ... ---- | df_i | dx_1 dx_2 dx_N | jac = ---- = | | dx_j | . . . . | | . . . . | | . . . . | | | | df_N df_N df_N | | ---- ---- ... ---- | | dx_1 dx_2 dx_N | The fourth argument is optional, and may be used to specify a set of times that the ODE solver should not integrate past. It is useful for avoiding difficulties with singularities and points where there is a discontinuity in the derivative. After a successful computation, the value of ISTATE will be 2 (consistent with the Fortran version of LSODE). If the computation is not successful, ISTATE will be something other than 2 and MSG will contain additional information. You can use the function `lsode_options' to set optional parameters for `lsode'. Examples =========================================================================== 1. Solve the following system on [0,5] using Octave. x' = y * sin(t) y' = x + cos(t) x(0) = 1 y(0) = 0 >> f = @(x,t) [x(2)*sin(t), x(1)*cos(t)]'; >> x0 = [1, 0]'; >> t = linspace(0, 5, 501); >> x = lsode(f, x0, t); >> plot(t, x(:,1), 'r-', t, x(:,2), 'b-' ) >> legend('x', 'y', 'location', 'southwest') >> xlabel('t') >> title('Solution of System') 2. Solve y''' + 2y'' - y' - 2y = exp(t) 0<=t<=3 subject to y(0) = 1, y'(0) = 2, y''(0) = 0. >> f = @(x,t) [x(2), x(3), exp(t) + dot([2,1,-2],x)]; >> t = linspace(0, 3, 31)'; >> x0 = [1,2,0]'; >> x = lsode(f, x0, t); >> plot(t, x(:,1)) >> xlabel('t') >> ylabel('y') >> title('Solution')