Info
Questa domanda è chiusa. Riaprila per modificarla o per rispondere.
Ode45 problem with equations
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
John D'Errico
il 19 Mag 2019
Modificato: John D'Errico
il 19 Mag 2019
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.
Risposte (1)
Star Strider
il 19 Mag 2019
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 Commenti
Star Strider
il 19 Mag 2019
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.
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!