Unable to solve stiff differential equation

Question

AJMIT KUMAR 2023-4-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1948648-unable-to-solve-stiff-differential-equation

评论： Sam Chak 2023-4-19

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更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2023-4-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1948648-unable-to-solve-stiff-differential-equation#answer_1217518

在 MATLAB Online 中打开

Hi @AJMIT KUMAR,

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.

% 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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sam Chak 2023-4-19

在 MATLAB Online 中打开

Hi @AJMIT KUMAR

Just to add something; your nonlinear system is unstable

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

.

If the following example, the sign of N is switched, and

stays within

.

Note: I didn't perform rigorous mathematical proof to show that the nonlinear system is stable for

. It is just my intuitionistic logic.

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))

ans = 5.2387

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

请先登录，再进行评论。

Answer 2

Oguz Kaan Hancioglu 2023-4-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1948648-unable-to-solve-stiff-differential-equation#answer_1217533

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.

https://www.mathworks.com/matlabcentral/answers/92961-how-do-i-use-a-fixed-step-size-with-ode23-and-ode45-in-matlab