In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
Lslack = 0.223
Umax = 0.04
Lceopt = 0.093
width  = 0.63*Lceopt
Fmax = 7400
a = 0.25*Fmax
b = 0.25*10*Lceopt

Lce = 0.087
dt = 0.001
time = np.arange(0,2.99,dt)
data = np.zeros((len(time),2))
In [3]:
for i, t in enumerate(time):

    if (t<=1):
        Lm = 0.31 
    elif (t>1 and t<2):
        Lm = 0.31 - 0.04*(t-1)
    else:
        Lm = 0.27
        
    Lsee = Lm - Lce
    
    if (Lsee < Lslack):
        F = 0
    else:
        F = Fmax*((Lsee-Lslack)/(Umax*Lslack))**2
        
       
    if (Lce < Lceopt):
        Fkpe = 0
    else:
        Fkpe = Fmax*((Lce-Lceopt)/(Umax*Lceopt))**2
        
    F0 = max([0, Fmax*(1-((Lce-Lceopt)/width)**2)])

    if (F > F0):
            print('Erro')
    Lcedot = -b*(F0-F+Fkpe)/(F+a)
    Lce = Lce + dt*Lcedot

    data[i,0] = t
    data[i,1] = F
In [4]:
plt.plot(data[:,0], data[:,1])
Out[4]:
[<matplotlib.lines.Line2D at 0x7d347b8>]
In [5]:
Lce
Out[5]:
0.042481884299063612
In [6]:
Lm
Out[6]:
0.27
In [ ]: