Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How Define ODE With Piecewise Function? #52

Open
FFFxueGawaine opened this issue Jun 28, 2024 · 0 comments
Open

How Define ODE With Piecewise Function? #52

FFFxueGawaine opened this issue Jun 28, 2024 · 0 comments

Comments

@FFFxueGawaine
Copy link

FFFxueGawaine commented Jun 28, 2024

Excuse me, I want to solve a single DOF oscillator with a nonlinear spring, whose nonlinear stiffness is modelled as Piecewise Function. How I define ODE?

I try use sympy.Piecewise. While there were no errors and warnings reported, the results appeared to be wrong.

I appended a picture of the correct result calculated by scipy.solve_ivp. The ODE is

p = [1, 0.1, 100, 1000,  0]
def impact_1dof(t, x, p):
    return [x[1],
            1 / p[0] * ( - p[1] * x[1] - p[2] * x[0] - p[3] * (x[0] > p[4]) * (x[0] - p[4]))]

The reslut is correct as follow.

微信图片_20240628184310

Next is code using sunode, and its reslut.

I need your help, if you understand what's going on. THANKS.

 
def SDOF(t, y, p):
 
   
    # nonlinear_damper = pt.switch(y.x > p.delta1, p.cnl * y.v, 0.0)
    return {
        'x':  y.v,
        'v': -1/p.m*(p.k*y.x+p.c*y.v +p.knl* sp.Piecewise((0, y.x <= 0), (y.x, y.x > 0)))
    }


problem = SympyProblem(
    params={
        'm': (),
        'k': (),
        'c': (),
        'knl':(),
      
       
           },
    states={
        'x': (),
        'v': (),
    },
    rhs_sympy=SDOF,
    derivative_params=[
        ('m',),
        ('c',),
        ('k',),
        ('knl',),
         
    ],
)

solver = Solver(problem, sens_mode=None, solver='BDF')
tvals = np.linspace(0, 5,10000)
 
y0 = np.zeros((), dtype=problem.state_dtype)
y0['x'] = -0.1
y0['v'] = 0
# We can also specify the parameters by name:
solver.set_params_dict({
    'm': 1 ,
    'k': 100,
    'c': 0.1,
    'knl': 100,
    
})
output = solver.make_output_buffers(tvals)
solver.solve(t0=0, tvals=tvals, y0=y0, y_out=output)

# We can convert the solution to an xarray Dataset
solx= solver.as_xarray(tvals, output).solution_x 
solv = solver.as_xarray(tvals, output).solution_v
plt.figure()
plt.plot(solx,solv)

微信图片_20240628184830

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant