Info
此问题已关闭。 请重新打开它进行编辑或回答。
Ode45 problem with equations
1 次查看(过去 30 天)
显示 更早的评论
I'm having problems with properly typing in four 2nd order differential equations into ode.45 function.
The code looks like this:
load parameters
t0=0:0.01:10;
X01=zeros(8,1);
for k = 1:length(xr)
[t1,X2] = ode45(@(t,x)f1(t,x,dxr(k),dyr(k),xr(k),yr(k)),t0,X01);
end
And the function
function dX=f1(t,x,dxr,dyr,xr,yr)
m1=10;
mn=10;
p = 1;
R = 0.01;
u = 0.7;
ks = 1.5 * 10^8;
g=9.81;
F1r = ((yr - x(6))^p)*ks*(1 - ((1 - (R^2))/2)*(1 - sign(yr - x(6))*sign(dyr - x(2))));
T1r = -u * F1r * sign(x(1) - dxr);
F21 = ((x(6) - x(8))^p)*ks*(1 - ((1 - R^2)/2)*(1 - sign(x(6) - x(8))*sign(x(2) - x(4))));
T21 = -u * F21 * sign(x(3) - x(1));
M=zeros(8,8);
M(1,1) = m1;
M(2,2) = m1;
M(3,3) = m1;
M(4,4) = m1;
M(5,5) = 1;
M(6,6) = 1;
M(7,7) = 1;
M(8,8) = 1;
dx1=x(1,1);
dy1=x(2,1);
dx2=x(3,1);
dy2=x(4,1);
x1=x(5,1);
y1=x(6,1);
x2=x(7,1);
y2=x(8,1);
Q = [T1r - T21; F1r-m1*g-F21; T21; F21-m1*g; dx1; dy1; dx2; dy2];
dX=inv(M)*Q;
end
It just gets stuck and doesn't execute single loop on ode45.
2 个评论
John D'Errico
2019-5-19
编辑:John D'Errico
2019-5-19
I itried looking at your code. The very first thing I would do is to execute your function F1, at the initial values. LOOK AT WHAT IT RETURNS. Does it execute at all? What happens? But that fails of course, because we don't have those other parameters, thus
%dxr, dyr, xr, yr are defined earlier parameters
That makes it impossible to execute your code to test it out.
If you want help, then make it easy for someone to help you.
回答(1 个)
Star Strider
2019-5-19
If your ‘f1’ function executes, then at least it is working.
I get the impression that it is ‘stiff’, particularly because of the value of ‘ks’.
The solution is likely to be to use ode15s instead of ode45.
2 个评论
Star Strider
2019-5-19
The warning is telling you that there is a singularity (infinite result) at about that time. Assuming your lower time limit is 0, set:
t0 = [0 0.001];
to see what your ‘f1’ function is doing up to that time.
Finding out the cause of the singularity and correcting it is something only you can do.
此问题已关闭。
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!