清空筛选条件
清空筛选条件

Problem when solving ODEs with ode45

1 次查看（过去 30 天）
Christo van Rensburg
评论： Ameer Hamza 2020-4-28
Good day,
I have the following code where I simultaneously solve three second-order differential equations defined in 'ode3'.
My problem is that in the second last line in 'ode3', the denominator is Y3(3). This is the variable 'r' and is in the direction of a spring and stretches as a mass swings around at its endpoint. When I run the script I get NaN in my Y3 matrix, which then leads me to believe I'm dividing by zero some time. What I don't understand is why would that value become zero?
Any help is appreciated!
m1 = 4; %kg
m2 = 6; %kg
L = 1.5; %m
k = 100; %N/m
g = 9.81; %m/s2
% --------------------------------------------------------------
% Force parameters
% --------------------------------------------------------------
Fo = 100; %N
tf = 1;
F = @(t) (t<=tf)*Fo*sin(2*pi*t/tf);
ic3 = [0, 0, 0, 0, L, 0];
ode3 = @(t, Y3) [Y3(4);
Y3(5);
Y3(6);
(F(t) - m2*k*(L-Y3(3))*sin(Y3(2)))/m1;
(-F(t)*cos(Y3(2)) + m2*k*(L-Y3(3))*sin(Y3(2))*cos(Y3(2)) - 2*m1*Y3(6)*Y3(5) - m1*g*sin(Y3(2)))/(m1*Y3(3));
(-F(t)*sin(Y3(2)) + m2*k*(L-Y3(3))*sin(Y3(2))^2 + m1*Y3(3)*Y3(5)^2 + m1*g*cos(Y3(2)) + m1*k*(L-Y3(3)))/m1];
[t3, Y3] = ode45(ode3, [0, 3], ic3);

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

采纳的回答

Ameer Hamza
Ameer Hamza 2020-4-28
In your initial condition, you have
ic3 = [0, 0, 0, 0, L, 0];
which implies that at initial point, Y(3) = 0. But in this equation
(-F(t)*cos(Y3(2)) + m2*k*(L-Y3(3))*sin(Y3(2))*cos(Y3(2)) - 2*m1*Y3(6)*Y3(5) - m1*g*sin(Y3(2)))/(m1*Y3(3))
m1*Y3(3) comes in the denominator. At the initial point, you get a 0/0 situation. If you try some other initial condition, say
ic3 = [0, 0, 1, 0, L, 0];
You will get a finite answer.
You will need to check you system of ODE, if it is correctly written, and which value of initial conditions are allowed.
  2 个评论
Christo van Rensburg
Christo van Rensburg 2020-4-28
Cool thanks so much!
It was a silly mistake from my side where the initial conditions where incorrectly allocated.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by