Ode45 problem with equations

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.

[John D'Errico]

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.

[buszcini]

That is my bad sir, I have updated the code and I put these parameters as attachment

Answer 1

[Star Strider]

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.

[buszcini]

I've tried ode15s but there is this warning

Warning: Failure at t=1.169317e-03. Unable to meet integration tolerances without reducing the step size below the smallest value

allowed (3.469447e-18) at time t.

[Star Strider]

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.