solve_ivp
Fixedstep IVP solvers for solving vectorvalued initial value problems.
Back to IVP Solver Toolbox Contents.
Contents
Syntax
[t,y] = solve_ivp(f,[t0,tf],y0,h) [t,y] = solve_ivp(f,{t0,E},y0,h) [t,y] = solve_ivp(__,method) [t,y] = solve_ivp(__,method,wb)
Description
[t,y] = solve_ivp(f,[t0,tf],y0,h) solves the IVP defined by f(t,y) from t0 until tf using the classic RungeKutta fourthorder method with an initial condition y0 and step size h.
[t,y] = solve_ivp(f,{t0,E},y0,h) does the same as the syntax above, but instead of terminating at a final time tf, the solver terminates once the event function E(t,y) is no longer satisfied.
[t,y] = solve_ivp(...,method) can be used with either of the syntaxes above to specify the integration method (i.e. the function will use the specified integration method, instead of the classic RungeKutta fourthorder method used by the previous two syntaxes). Options for the integration method are listed in the "Input/Output Parameters" section.
[t,y] = solve_ivp(...,method,wb) can be used with any of the syntaxes above to define a waitbar. If wb is input as true, then a waitbar is displayed with the default message 'Solving IVP...'. To specify a custom waitbar message, input wb as a char array storing the desired message.
Input/Output Parameters
Variable  Symbol  Description  Format  
Input  f  multivariate, vectorvalued function () defining the ordinary differential equation
 inputs to f are the current time (t, 1×1 double) and the current state vector (y, p×1 double)  output of f is the state vector derivative (dydt, p×1 double) at the current time/state 
1×1 function_handle 

t0  initial time  1×1 double 

tf  final time  1×1 double 

C  event function ()
 inputs are the current time (t, 1×1 double) and the current state vector (y, p×1 double)  output is a 1×1 logical (true if solver should continue running, false if solver should terminate) 
1×1 function_handle 

y0  initial condition  p×1 double 

h  step size  1×1 double 

method    (OPTIONAL) integration method (defaults to 'RK4')
 'Euler' (Euler 1storder method)  'RK2' (midpoint method)  'RK2 Heun' (Heun's 2ndorder method)  'RK2 Ralston' (Ralston's 2ndorder method)  'RK3' (Kutta's 3rdorder method)  'RK3 Heun' (Heun's 3rdorder method)  'RK3 Ralston' (Ralston's 3rdorder method)  'SSPRK3' (strong stability preserving 3rdorder method)  'RK4' (classic RungeKutta 4thorder method)  'RK4 Ralston' (Ralston's 4thorder method)  'RK4 3/8' (3/8rule 4thorder method)  'AB2' (AdamsBashforth 2ndorder method)  'AB3' (AdamsBashforth 3rdorder method)  'AB4' (AdamsBashforth 4thorder method)  'AB5' (AdamsBashforth 5thorder method)  'AB6' (AdamsBashforth 6thorder method)  'AB7' (AdamsBashforth 7thorder method)  'AB8' (AdamsBashforth 8thorder method)  'ABM2' (AdamsBashforthMoulton 2ndorder method)  'ABM3' (AdamsBashforthMoulton 3rdorder method)  'ABM4' (AdamsBashforthMoulton 4thorder method)  'ABM5' (AdamsBashforthMoulton 5thorder method)  'ABM6' (AdamsBashforthMoulton 6thorder method)  'ABM7' (AdamsBashforthMoulton 7thorder method)  'ABM8' (AdamsBashforthMoulton 8thorder method) 
char  
wb    (OPTIONAL) waitbar parameters (defaults to false)
 input as "true" if you want waitbar with default message displayed  input as a char array storing a message if you want a custom message displayed on the waitbar 
char or 1×1 logical 

Output  t  time vector  (N+1)×1 double 

y  solution matrix  (N+1)×p double 
Time vector and solution matrix:
Note
 The nth row of stores the transpose of the state vector (i.e. the solution) corresponding to the nth time in . This convention is chosen to match the convention used by MATLAB's ODE suite.
Example #1: Time detection.
Consider the Lorenz system, with , , and :
Plot the solution of this system for in the interval .
First, we can rewrite this system in vector form by letting , , , and .
Defining the system and its initial condition in MATLAB,
% Lorenz parameters rho = 28; sigma = 10; beta = 8/3; % Lorenz equations in vector form f = @(t,x) [sigma*(x(2)x(1)); x(1)*(rhox(3))x(2); x(1)*x(2)beta*x(3)]; % initial condition x0 = [0; 1; 1.05];
Solving the system for in the interval using the classic RK4 method with a time step (i.e. step size) of ,
[t,x] = solve_ivp(f,[0,100],x0,0.001);
Plotting the solution,
figure; plot3(x(:,1),x(:,2),x(:,3)); view(45,20); grid on; xlabel('$x$','Interpreter','latex','FontSize',18); ylabel('$y$','Interpreter','latex','FontSize',18); zlabel('$z$','Interpreter','latex','FontSize',18);
Example #2: Event detection.
Consider the initial value problem
Find the solution until .
First, let's define our ODE () and initial condition in MATLAB.
f = @(t,y) y; y2 = 3;
Next, let's define the event function, . Since we want to continue solving until , we want the solver to run while ; this forms our condition. Therefore,
E = @(t,y) y <= 10;
Solving for using the Euler method with a step size of ,
[t,y] = solve_ivp(f,{2,E},y2,0.01,'Euler');
Plotting the solution,
figure; plot(t,y,'LineWidth',1.5); grid on; xlabel('$t$','Interpreter','latex','FontSize',18); ylabel('$y$','Interpreter','latex','FontSize',18);
Example #3: Backward integration (time detection case).
Consider the same ODE as in Example #2, but we are now given its value at .
Find . Then, confirm your result by solving the same ODE from to , using as the initial condition.
In this case, we know at , and want to find at a previous time . To do this, we need to solve the ODE backwards, from to . First, let's define our ODE () and initial condition in MATLAB.
f = @(t,y) y; y20 = 50;
Solving for from to using the ABM8 method with a step size of ,
[t,y] = solve_ivp(f,[20,10],y20,0.001,'ABM8');
The solution for corresponding to will be located at the last element of the solution matrix, since is stored in the last element of the time vector. Therefore, is
y10 = y(end)
y10 = 0.002269996488128
Confirming our result by solving the same ODE but from to and using our result for as the initial condition,
[t,y] = solve_ivp(f,[10,20],y10,0.001,'ABM8');
y20 = y(end)
y20 = 49.999999999999226
Example #4: Backward integration (eventdetection case).
Once again, consider
Find using an IVP solver with a event function (i.e. eventdetection).
First, let's define our ODE () and initial condition in MATLAB.
f = @(t,y) y; y20 = 50;
The event that terminates the solver is when . Therefore, we define the event function as
E = @(t,y) t > 10;
In the timedetection case (Example #3) we input [t0,tf] = [20,10], so the solver knew to integrate backwards in time since t0 > tf. Consequently, internally, the solver made the step size negative. However, for the eventdetection case, just given t and E(t,y), the solver won't know to use a negative step size to integrate backwards in time. Therefore, we must manually specify a negative step size. Solving for ,
[t,y] = solve_ivp(f,{20,E},y20,0.001,'ABM8');
y10 = y(end)
y10 = 0.002272267619989
Note that this is the same result we obtained earlier in Example #3.