Issue with integration tolerances

조회 수: 12 (최근 30일)
Michael
Michael 2025년 2월 25일
댓글: Walter Roberson 2025년 2월 25일
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
Torsten
Torsten 2025년 2월 25일
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 2025년 2월 25일
I tried it symbolically, taking the given equations and using dsolve()... it said that it could not find the answer.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by