Adaptive step size solver for an ODE

This question involves both mathematics and MATLAB. I am trying solve an ordinary differential equation numerically, dy/dt=10*exp(-(t−2)^2/(2*(0.075)^2)−0.6y with the initial value y(0)=0.5 and initial step size h=0.3 between t=0 and t=4. In my code I implement an embedded Runge-Kutta formulation: Calculating two tentative estimates using RK2 and RK3 formulae with the same step size and the difference between them is the error for RK2. Then I compare this difference (will simply called 'error') with a value (an acceptable error) which actually is not a set number but is obtained through 10^-3∗|ycurrent| or 10^−6 whichever is larger (ycurrent is the present estimate) If the error is greater than the acceptable error I update step size with hnew=hold(acceptableerrorerror)13 and if the error is less than or equal to the acceptable error, the latest step size holdis used to estimate the function value and the independent variable t is increased by that amount.

t(1)=0; y(1)=0.5; h(1)=0.3;
i=1; j=1;
while t<=4

 k1=f(t(i),y(i));
 k2=f(t(i)+0.5*h(j),y(i)+0.5*h(j)*k1);
 k3=f(t(i)+0.75*h(j),y(i)+0.75*h(j)*k2);    % finding k's For RK3 formula

 y(i+1)=y(i)+(h(j)/9)*(2*k1+3*k2+4*k3);     % estimate using RK3

 t(i+1)=t(i)+h(j);                          % calculating the next t to find k4
 k4=f(t(i+1),y(i+1));
 E(i+1)=abs((h(j)/72)*(-5*k1+6*k2+8*k3-9*k4));   % local truncation error, difference between RK2 and RK3 estimations
 y(i+1)=y(i)+(h(j)/9)*(2*k1+3*k2+4*k3)-(h(j)/72)*(-5*k1+6*k2+8*k3-9*k4);
 epsilon_s=max(10^-3*abs(y(i)),10^-6);
    if  E(i+1)>epsilon_s
        h(j+1)=0.8*h(j)*(epsilon_s/E(i+1))^(1/3);   % updating step size
        t(i+1)=t(i)+h(j);
        j=j+1;
    else
        t(i+1)=t(i)+h(j);                           % if actual error falls below or is equal to epsilon_s, current step size is accepted
        h(j+1)=h(j)/(0.8*(epsilon_s/E(i+1))^(1/3)); % I need to increase step size here
        i=i+1;
        j=j+1;
        t(i+1)=t(i)+h(j);                           % this line may make no sense but i need the next t for k4 above
    end
end

The local truncation error using RK2 is the absolute value of the difference between the estimations made by RK2 and RK3 formulae. When the error criterion is satisfied I proceed with the y which is estimated by RK2. The line

y(i+1)=y(i)+(h(j)/9)*(2*k1+3*k2+4*k3)-(h(j)/72)*(-5*k1+6*k2+8*k3-9*k4);

is RK2 estimate of y(i+1). It is actually the difference between RK3 estimate and the local truncation error for using RK3.

I have two issues here: The script runs indefinitely and h (step size) quickly goes to 0. The second is that I do not feel confident with the line

h(j+1)=h(j)/(0.8*(epsilon_s/E(i+1))^(1/3)); % I need to increase step size here

because here I use the inverse of the

h(j+1)=0.8*h(j)*(epsilon_s/E(i+1))^(1/3)

formula to increase step size. (Step size is made larger here beacuse the error requirement is just met there is room for me to go on with a larger step size)

I ask for maybe a different way of increasing my step size.

By the way, f is a function m-file

function z=f(t,y);
z=10*exp((-(t-2)^2)/(2*0.075^2))-0.6*y;

[Jan]

This is the standard to control the stepsize:

h(j+1) = 0.8*h(j)*(epsilon_s/E(i+1))^(1/3);

In case of rejected steps the stepsize is reduced by 50% usually:

h(j+1) = h(j) / 2;

This is applied multiple times:

t(i+1)=t(i)+h(j);

Remember, that 10^-3 is an expensive power operation, while 1e-3 is a cheap constant.

[Ali Kiral]

Jan,can I use the same formula

h(j+1) = 0.8*h(j)*(epsilon_s/E(i+1))^(1/3)

for both increasing and decreasing the step size? or does h(j+1) always happen to be smaller than h(j)?

and

k4=f(t(i+1),y(i+1))

is the reason of multiple application of

t(i+1)=t(i)+h(j)

Do you suggest I use 1e-3 and 1e-6 instead of powers?

Answer 1

[Sam Chak]

Hi @Ali Kiral

After fixing the code, you can compare your result with the result by ode45().

A little analysis shows that the system

is stable because the Gaussian disturbance

dissipates over time , making the system behaves like an exponential decay

.

% Solving the ODE using ode45
odefcn = @f;
tspan  = linspace(0, 10, 1001);
y0     = 0.5;                           % initial value
[t, y] = ode45(odefcn, tspan, y0);

% disturbance 
d = 10*exp((-(t - 2).^2)/(2*0.075^2));  

% Viewing the effect of disturbance d(t) on the system response y(t)
plot(t, [d y]), grid on, xlabel('t')
legend('Disturbance, d(t)', 'Sys Response, y(t)')

Also, the analytical solution exists and it can expressed in terms of the Gauss error function erf(), which is a non-elementary integral.

% Analytical Solution of the System
syms y(t) 
eqn     = diff(y,t) == - 0.6*y + 10*exp((-(t - 2)^2)/(2*0.075^2));
cond    = y(0)      ==   0.5;
ySol(t) = dsolve(eqn, cond)

ySol(t) = 

% Differential equation
function dydt = f(t, y);
    dydt = 10*exp((-(t - 2)^2)/(2*0.075^2)) - 0.6*y;
end