How can I use an IF statement to define a variable when solving a Differential equation

Dear All,

I'm trying to solve a simple oscillating mass spring system and I would like to have the spring stiffness depending on the position of the mass.

Here I define the Equation of motion:

function F= MSI(t,x,flag,m,c,g)

F = [x(2); -(c/m)*x(1)- g];

And the in an other file I let a solver solve it:

m=1;
g=9.81;

if x(:,1)>0
    c=0;
else
    c=10000;
end

tf=15;
ts=0.01;
time=(0:ts:tf);
N=length(time);

y0=1;
yd0=0;
initial= [y0 yd0];

tic;
options=odeset('AbsTol',1.e-24,'RelTol',1.e-12);  
[t x]=ode113('MSI',time,initial,options,m,c,g);
elapsed_time_sec=toc;

f1=figure;
 plot(t(:,1),x(:,1))
     ylabel('y','FontSize',11)

You can see that I try to let spring stiffness c depend on position y, or in this case x(:,1). This however doesn't work.

How can I feedback y which is written down in x(:,1) to determine c and feed it back into the EOM and solver?

Best regards,

Gerben Smit

[Jan]

What exactly does "it doesn't work" mean? Do you get an error message, a warning or do the results differ from your expectations? It is much easier to solve a problem, when we do not have to guess it at first.

Answer 1

[Jan]

An absolute tolerance of 1e-24 is ridiculously small. Even a relative tolerance of 1e-12 can hurt the integrator seriously. Much greater tolerances are strongly recommended for a successful integration.