![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1359043/image.png)
Unable to solve stiff differential equation
3 次查看(过去 30 天)
显示 更早的评论
I have solved many nonlinear differential using ode45 previously.
For the current stiff nonlinear differential equation, I used ode45, ode15s, ode23s, ode23t and ode23tb but unable to get final results.
I got constant warnings which is mentioned below.
Warning: Failure at t=8.048300e+00. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (2.859330e-14) at time t.
Warning: RelTol has been increased to 2.22045e-14.
The code is below.
clear all;clc;
op=odeset('relTol',1e-14,'absTol',1e-14);
global psi D B H M epsilon N
epsilon=0.01;
M=epsilon*1;
H=epsilon*1;
psi=5;
D=epsilon*1;
B=epsilon*0.5;
N=epsilon*0.17;
w= 0.1:0.01:5;
for i=1:1:length(w)
THETA=w(i);
[t,y]=ode23s(@(t,x)F(t,x,THETA),[0 20],[0 0],op);
end
function dy = F(t,y,THETA)
global psi D B H M N
dy(1) = y(2);
dy(2) = psi.*cos(THETA.*t)-(1-M.*sin(THETA.*t)).*y(1)-B.*y(2)-(H.*(cos(THETA.*t))).*((y(1)).^2)-(D.*(sin(THETA.*t))-N)*((y(1))^3);
dy=dy';
end
0 个评论
回答(2 个)
Sam Chak
2023-4-17
It appears that the system is inherently unstable for
. If you simulate the system for a longer time, then the system will explode. Small value of Θ takes longer to explode.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1359043/image.png)
% op = odeset('relTol', 1e-14, 'absTol', 1e-14);
global psi D B H M epsilon N
epsilon = 0.01;
M = epsilon*1;
H = epsilon*1;
psi = 5;
D = epsilon*1;
B = epsilon*0.5;
N = epsilon*0.17;
w = 0.19:0.001:0.195; % theta 0.198 will stop integrating at around 20 sec
% w = zeros(1, 2); % stable
% w = 0.01:0.001:0.014; % stable
for i = 1:1:length(w)
THETA = w(i);
[t, y] = ode45(@(t,x) F(t, x, THETA), [0 20], [0 0]);
plot(t, y), hold on
end
hold off, grid on, xlabel('t, (sec)')
function dy = F(t, y, THETA)
global psi D B H M N
dy(1) = y(2);
dy(2) = psi*cos(THETA*t) - (1 - M*sin(THETA*t))*y(1) - B*y(2) - H*cos(THETA*t)*y(1)^2 - (D*sin(THETA*t) - N)*y(1)^3;
dy = dy';
end
1 个评论
Sam Chak
2023-4-19
Hi @AJMIT KUMAR
Just to add something; your nonlinear system is unstable
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1360993/image.png)
because the coefficient of the cubic growth term
is "POSITIVE". In other words, if y gets large, it rapidly "adds" energy to
. On the contrary, if you make
, then the system will be stable because it rapidly "dissipates" energy in
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1360998/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361003/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361008/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361013/image.png)
If the following example, the sign of N is switched, and
stays within
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361018/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361023/image.png)
Note: I didn't perform rigorous mathematical proof to show that the nonlinear system is stable for
. It is just my intuitionistic logic.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361028/image.png)
global psi D B H M epsilon N
epsilon = 0.01;
M = epsilon*1;
H = epsilon*1;
psi = 5;
D = epsilon*1;
B = epsilon*0.5;
N = - epsilon*0.17; % <--- this coefficient must be negative
w = 0.19:0.001:0.195; % theta 0.198 will stop integrating at around 20 sec
% w = zeros(1, 2); % stable
% w = 0.01:0.001:0.014; % stable
for i = 1:1:length(w)
THETA = w(i);
[t, y] = ode45(@(t,x) F(t, x, THETA), [0 3000], [0 0]);
plot(t, y(:,1)), hold on
end
hold off, grid on, xlabel('t, (sec)'), ylabel('y(t)')
title({'Plot of the solution for $y_{1}(t)$'}, 'interpreter', 'latex', 'fontsize', 16)
idx = find(t>2000);
ymax = y(:,1);
max(ymax(idx))
function dy = F(t, y, THETA)
global psi D B H M N
dy(1) = y(2);
dy(2) = psi*cos(THETA*t) - (1 - M*sin(THETA*t))*y(1) - B*y(2) - H*cos(THETA*t)*y(1)^2 - (D*sin(THETA*t) - N)*y(1)^3;
dy = dy';
end
Oguz Kaan Hancioglu
2023-4-17
I think matlab solver couldn't solve the nonlinear probrem within your tolerances. You can use fixed step size which is very small dt values (e-5, or e-6) to solve your problem.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!