How to fix the time step in ODE45

I want to fix the time step for my ode45 function. The code which I am using is as follows:

    dt = 0.02;
    tf = 600;
    t = dt:dt:tf;
    y0 = zeros(14,1);
    [tout,yout] = ode45(@OC3_odefull,t,y0);

When I run my code, I have no control over the time step size and ode45 uses an adaptive time step. Is there any way that I can force ode45 to use the time step that I want?

[Jan]

The "45" means, that each step is calculated with an order 4 and order 5 method. The difference between the results is used to control the step size. It is not meaningful to run an ODE45 method with a fixed step size.

[Abla Nevame Edzenunye]

Thank you!!!

Answer 1

[Steven Lord]

Why do you want to fix the time step?

If you want to find the solution of the system of ODEs at specific times, you don't need to control the time step to do that. Specify the time span vector as a vector with more than two elements and the ODE solver will return the solution at the specified times. [Note that this affects the times at which the solver returns the solution, not the times used by the solver internally.]

If you have a differential equation where the value of the right-hand side depends upon the value of the solution at earlier times (and you're trying to ensure the solver computes the solution at those earlier times) you don't want to use the ordinary differential equation solvers. Use the delay differential equation solvers instead, like dde23 or ddesd.

If there's a different reason you want to control the time step, please say more about that reason.

[Vahid Jahangiri]

I am trying to load a file inside my code which is based on the time span I want.

[Walter Roberson]

A typical way of handling that situation is to interpolate the values from the file according to the current time.

[Jan]

But for interpolating keep in mind, that Matlab's ODE integrators require a smoothly differntiable function. So do not use a linear interpolation, because it will confuse the step size controller. A cubic interpolation is fine.

[Walter Roberson]

Note that it is not advisable to load the file each time. Instead, load it before the ode*() call, and pass the data in as parameters.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Vishesh Mangla]

Sorry, it doesn't seem to work. I tried using tspan=linspace(.1, 1, 10) and I got 56 values instead of 10.

ans =

Columns 1 through 14

0.1000 0.1167 0.1335 0.1664 0.2088 0.2304 0.2521 0.2830 0.3126 0.3273 0.3421 0.3590 0.3751 0.3888

.......

Columns 43 through 56

0.8094 0.8245 0.8422 0.8553 0.8684 0.8839 0.9009 0.9165 0.9296 0.9419 0.9561 0.9728 0.9896 1.0000

I can see close enough values amongst these "times" but I want to plot this data and it is quite cumbersme to do all those slicing operations to select whatever's required.

[Steven Lord]

Can you show your ode45 call (including the code with which you define all the variables used in that call)?

[Walter Roberson]

That ans that is shown, is it the output from linspace() ? Or is it the output times from an ode45 call?

The entries are not equally spaced or nearly so; linspace() would never produce that output.

[Vishesh Mangla]

v 2018a

This is the whole code and you can replicate by just copying. Just try obj.temperature function which takes argument as a vector, so try linspace or [0 .1 .2 ] but still instead of 3 points on the time axis there is all that stuff making the plot dirty.

FiniteElement.m a class

classdef FiniteElement
    %FINITEELEMENT Summary of this class goes here
    
    properties
        simpson_nodes 
        N%number of nodes
        H% step in x direction
        domain % 0 to 1
        beta% hyperparameter
        
        %initial conditions
        alpha
        
        %analytic solution
        exact
        
        f%f(x, t)
        
        
        %Define boundary conditions
        g0  % alpha0
        g1  %alphaN
        g0dash  %alpha'0
        g1dash  %alpha'N
        
        %stiffness matrices
        A
        D
        % time to calculate at
        time_at
    end
    
    methods
        
        function obj = FiniteElement(N, beta)
            
            % Constructor initialize all variables
            %Parameters Initialization
            obj.N = N;
            obj.beta = beta;
            obj.H = 1 / N; % step in x direction
            obj.domain = linspace(0, 1, N + 1);
            obj.f = @(x, t) (1 - pi^2 * obj.beta ) * exp(t) * sin(pi * x);
            
            %Set initial conditions here
            obj.alpha = sin(pi * obj.domain(2:end - 1));
            
            %Set exact solution here
            obj.exact = @(t) exp(t) * obj.alpha;
            obj.simpson_nodes = 3*20*obj.N; % number of nodes to be used for
            %Define boundary conditions
            obj.g0 = @(t) 0;
            obj.g1 = @(t) 0;
            obj.g0dash = @(t) 0;
            obj.g1dash = @(t) 0;
            obj = obj.matrices_req();
        end
        
        function output = phi(obj, i, x)
            % basis functions see live script for better representation
            if i == 0
                
                if x >= 0 && x <= obj.H
                    output = -(x - 1 * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H  &&  x <= obj.N * obj.H
                    output = (x - (obj.N - 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if x >= (i - 1) * obj.H && x <= i * obj.H
                    output = (x - (i - 1) * obj.H) / obj.H;
                elseif x >= i * obj.H && x <= (i + 1) * obj.H
                    output = -(x - (i + 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = phidash(obj, i, x)  % derivative
            
            %derivative of phi
            if i == 0
                
                if x >= 0 && x <= obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H && x <= obj.N * obj.H
                    output = 1 / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if x >= (i - 1) * obj.H && x <= i * obj.H
                    output = 1 / obj.H;
                elseif x >= i * obj.H && x <= (i + 1) * obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = const(obj, i, t)       
            % constant in the matrix equation
            output = obj.g0dash(t) * simpsons(@(x) obj.phi(i, x) * obj.phi(0, x), 0, 1, obj.simpson_nodes) + ...
                obj.g1dash(t) * simpsons(@(x) obj.phi(i, x) * obj.phi(obj.N, x), 0, 1, obj.simpson_nodes) ...
                -obj.beta * obj.g0(t) * simpsons(@(x) obj.phidash(i, x) * obj.phidash(0, x), 0, 1, obj.simpson_nodes) + ...
                -obj.beta * obj.g1(t) * simpsons(@(x) obj.phi(i, x) * obj.phi(obj.N, x), 0, 1, obj.simpson_nodes)...
                +obj.beta*obj.phi(i, 1)*obj.phidash(obj.N, 1)*obj.g1(t)-obj.beta*obj.phi(i, 0)*obj.phidash(0, 0)*obj.g0(t) ;
        end
        
        function obj = matrices_req(obj)
            % stiffness matrices
            obj.A = zeros(obj.N - 1);
            obj.D = zeros(obj.N - 1);
            
            for i = 1:obj.N - 1
                
                for j = 1:obj.N - 1
                    
                    if ismember(abs(i - j), [0, 1])
                        obj.A(i, j) = simpsons(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1, obj.simpson_nodes);
                        obj.D( i, j) = simpsons(@(x) obj.phidash(i, x) * obj.phidash(j, x), 0, 1, obj.simpson_nodes);
                    end
                    
                end
                
            end
            
        end
        
        
        function derivative_ = dydx(obj, t, y)
            % derivative used in ode45
            C = zeros(obj.N - 1, 1);
            F = zeros(obj.N - 1, 1);
            
            for i = 1:obj.N - 1
                C(i) = obj.const(i, t);
                F(i) = simpsons(@(x) obj.phi(i, x) * obj.f(x, t), 0, 1, obj.simpson_nodes);
                
            end
            
            derivative_ = obj.A \ (-(-obj.beta * obj.D * y + C) + F);
            
        end
        
        function output = exactsoln(obj, time)
            % give exactsolution the same structure as output from ode45
            output.x = time;
            output.y = [];
            
            for t =output.x
                output.y = [output.y obj.exact(t)'];
            end

            
        end
        function [ApproxTemp, ExactTemp] = temperature(obj, time_at)
            
            % main function to be called to get temperatures
            ApproxTemp = ode45(@(t, y) obj.dydx(t, y), time_at, obj.alpha);
            ExactTemp = obj.exactsoln(ApproxTemp.x);
        end 
    end
    
end

simpsons.m

function I = simpsons(f,a,b,n)

     h=(b-a)/n; 
     
     s = f(a) + f(b);
    for i=1:2:n-1  
        s = s+ 4*f(a+i*h);
        
    end
    for i=2:2:n-2  
        s = s + 2*f(a+i*h);
        
    end
    I = h/3 * s;
end

run.m

obj = FiniteElement(8, -1);
[A, E] = obj.temperature([0 .1 .2]); 
x = obj.domain;
y = A.x;
[xx, yy] = meshgrid(x, y);
subplot(2, 1, 1);
z = vertcat(arrayfun(obj.g0, A.x), A.y,arrayfun(obj.g1, A.x) )';
heatmap(x, y, z);
colormap(flipud(hot));
title("Approximate Solution");
xlabel("The rod as x-axis");
ylabel("Time");



subplot(2, 1, 2);
z = vertcat(arrayfun(obj.g0, E.x), E.y,arrayfun(obj.g1, E.x) )';
heatmap(x, y, z);
colormap(flipud(hot));
title("Exact Solution");
xlabel("The rod as x-axis");
ylabel("Time");

Run run.m

[Vishesh Mangla]

Sorry, mb. I just saw that there are couple of output formats and the one I was using was the structure output one. My apologies. I should have read the docs more thoroughly.

Answer 2

[Walter Roberson]

No, you cannot do that.

Officially, Mathworks only provides variable step solvers for MATLAB use.

Unofficially, see https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b#answer_107643 for the fixed step routines ode1 through ode5 .