Issue with integration tolerances

10 visualizaciones (últimos 30 días)
Michael
Michael el 25 de Feb. de 2025
Comentada: Walter Roberson el 25 de Feb. de 2025
I am trying to use this code to plot a series of Differential equations in Matlab, and it does give me plots however they are not spanning over the entire boundary of t that I am setting and I am getting this error
Warning: Failure at t=4.470044e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.
Can anyone explain why this is happening?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rudimentary script to solve an initial value problem consisting of the
% following system of two differential equations on the interval a<=t<=b
%
% x1'(t) = alpha*x1+beta*x2, x1(a) = x1_0
% x2'(t) = gamma*x1+delta*x2, x2(a) = x2_0
%
% Look for "*** Your entry i ***" where i=1,2,3,4,5,6 to make appropriate
% changes to solve your specific system
%
% Running the script as is outputs the unit circle as the trajectory with
% x1(t) = cos t and x2(t) = -sin t for the time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Housekeeping
clc ; % clear the command window
clear; % clear all variables
close all; % close all figures
% ODE45 options
options = odeset('AbsTol',1e-10,'RelTol',1e-6);
% Independent variable
a = 0; % left endpoint of t interval *** Your entry 1 ***
b = 20; % right endpoint of t interval *** Your entry 2 ***
h = 0.01; % stepsize *** Your entry 3 ***
tt = a:h:b;
% Initial value of dependent variables
x1_0 = 5; % initial value of x1 *** Your entry 4 ***
x2_0 = 1; % initial value of x2 *** Your entry 5 ***
x0 = [x1_0 x2_0];
% Solve the system
[t,soln] = ode45(@(t,x) rhs(t,x),tt,x0,options);
Warning: Failure at t=4.474016e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t.
% Plot solutions
figure(1); % time series of x1(t)
plot(t,soln(:,1));
xlabel('$t$','Interpreter','latex','FontSize',20);
ylabel('$x_1$','Interpreter','latex','FontSize',20,'Rotation',0)
grid on
figure(2); % time series of x2(t)
plot(t,soln(:,2));
xlabel('$t$','Interpreter','latex','FontSize',20);
ylabel('$x_2$','Interpreter','latex','FontSize',20,'Rotation',0);
grid on
figure(3); % phase portrait
plot(soln(:,1),soln(:,2));
xlabel('$x_1$','Interpreter','latex','FontSize',20);
ylabel('$x_2$','Interpreter','latex','FontSize',20,'Rotation',0)
grid on
axis square
% Right hand side of the system differential equations
function [xp] = rhs(t,x) % do not change this line
% *** Your entry 6 *** (start entering your system here)
alpha = 1.5;
beta = 1.1;
gamma = 2.5;
delta = 1.4;
x1prime = -alpha*x(1) + beta*x(1)*x(2);
x2prime = gamma*x(2)*(1-0.5*x(2)) + delta*x(1)*x(2);
% *** Your entry 6 *** (stop entering your system here)
xp = [x1prime ; x2prime]; % do not change this line
end % do not change this line
  2 comentarios
Torsten
Torsten el 25 de Feb. de 2025
Most probably, the solution to your ODE system has a singularity at t = 0.4474 ... (like tan(t) blows up at t = pi/2, e.g.) .
Walter Roberson
Walter Roberson el 25 de Feb. de 2025
I tried it symbolically, taking the given equations and using dsolve()... it said that it could not find the answer.

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Productos


Versión

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by