Problem with ODE solvers
2 views (last 30 days)
Show older comments
I want to numerically solve the following system of diff. eqns:
dx/dt = y
dy/dt = z
dz/dt = -y + x^2 - 0.025
where d is a constant.
The initial point is [0.5,0.5,0.5]
I have tried both stiff and nonstiff ODE solvers to no avail. The system always diverges to huge values and the following error pops up:
"Warning: Failure at t=2.660704e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t."
I have tried changing the eps value but it is not of much help.
A theoretical analysis shows that the system should not diverge. What do you guys suggest I do? The code that I used:
function [x,y,z] = system2(d ,initV, T, eps)
% d - system paramter
% INITV - initial point
% T - time interval
% EPS - ode solver precision
%
% Example.
% [X Y Z] = system2(1, [1.5 1.5 1.5], [0 25], 0.000001);
% plot3(X,Y,Z);
if nargin<2
eps = 0.000001;
T = [0 25];
initV = [1.5 1.55 1.5];
end
options = odeset('RelTol',eps,'AbsTol',[eps eps eps/100]);
[T,X] = ode23tb(@(T,X) F(T, X, d), T, initV, options);
plot3(X(:,1),X(:,2),X(:,3));
axis equal;
grid;
title('System2');
xlabel('X'); ylabel('Y'); zlabel('Z');
x = X(:,1);
y = X(:,2);
z = X(:,3);
return
end
function dx = F(T, X, d)
dx = zeros(3,1);
dx(1) = X(2);
dx(2) = X(3);
dx(3) = -X(2) + X(1)^2 - d;
return
end
0 Comments
Answers (1)
Amit
on 25 Jan 2014
What theoretical analysis you did which says that this does not diverges?
To me, this will diverge as the size of x goes bigger. Your initial value is > 1 for all X. There is no terms in all 3 differentials, that will keep it in bound. dX(3) will grow very fast (X(1)^2 will dominate very soon). This will make dX(2) grow bigger, making dx(1) bigger.
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!