# Adaptive step size solver for an ODE

28 views (last 30 days)
Ali Kiral on 2 Jan 2023
Answered: Sam Chak on 4 Jan 2023
This question involves both mathematics and MATLAB. I am trying solve an ordinary differential equation numerically, dy/dt=10*exp(-(t2)^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;
##### 2 CommentsShowHide 1 older comment
Ali Kiral on 2 Jan 2023
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?

Sam Chak on 4 Jan 2023
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

### Categories

Find more on Programming in Help Center and File Exchange

R2014a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!