Solving diff eq - piece wise

Python Help
1

im doing a physics lab and I’m trying to solve this system of equations:

(0.11551-0.002744Cv^2)/0.2025 = dv/dt
v=ds/dt
C=322/3300v when 3300v<10
C=10^(3.02-1.89*log(3300v)+0.411(log(3300v))^2-0.033(log(3300v))^3 when 10<3300v<20000
C=0.74 when 3300v>20000

v is velocity, C is drag coefficient, s is distance traveled, t is time
I’m trying to write a python code to find t at s=0.75. t is about 1.6, but i just want to get a more accurate answer. inital v=0, t=0, s=0. s is distance, v=velo, t=time. s, v, t, and C can’t be negative

here are some attempts (i’ve never coded in python before, i should probably learn some)

I’m having trouble incorporating C into the equation for dv/dt. I can calculate C in the 1st one, but I can’t implement it into the main equation, and in the 2nd one, I can’t implement a piece-wise C into dv/dt.

from Solving Second Order Differential Equations in Python - ePythonGuru

from matplotlib import pyplot as plt
from scipy.integrate import odeint
import numpy as np

def f(u, x):
    return [u[1], (0.11551 - 0.002744 * 0.74 * (u[1])**2) / 0.2025]

y0 = [0, 0]
xs = np.linspace(0, 10, 100)
us = odeint(f, y0, xs)
ys = us[:, 0]

plt.plot(xs, ys, '-')  # Plot the line
plt.plot(xs, ys, 'r*')  # Plot the points
plt.xlabel('x values')
plt.ylabel('y values')
plt.title('Updated Plot Title')
plt.show()

from chatgpt lol

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

# Function representing the modified differential equation
def model(v, t):
    # Calculate C based on the value of 3300v
    if 3300 * v < 10:
        C = 322 / (3300 * v)
    elif 10 < 3300 * v < 20000:
        C = 10**(3.02 - 1.89 * np.log(3300 * v) + 0.411 * (np.log(3300 * v))**2 - 0.033 * (np.log(3300 * v))**3)
    else:
        C = 0.74  
    
    # Your differential equation logic using v
    # For now, returning v itself without using C in the equation
    return 0.11551 / 0.2025

# Initial condition: v(0) = v0
v0 = 0.0  # Initial velocity value, change as needed

# Time points at which to solve for v (can be adjusted based on your requirements)
t_points = np.linspace(0, 10, 1000)  # Adjust the time range as needed

# Solve the modified differential equation
v_values = odeint(model, v0, t_points)

# Calculate s values using cumulative sum of v values and corresponding time points
s_values = np.cumsum(v_values.flatten()) * (t_points[1] - t_points[0])  # Integration sum

# Calculate C values based on v values
C_values = []
for v in v_values:
    v = v[0]
    if 3300 * v < 10:
        C = 322 / (3300 * v)
    elif 10 < 3300 * v < 20000:
        C = 10**(3.02 - 1.89 * np.log(3300 * v) + 0.411 * (np.log(3300 * v))**2 - 0.033 * (np.log(3300 * v))**3)
    else:
        C = 0  # Placeholder if conditions are not met
    C_values.append(C)

# Plotting t vs C, t vs s, and t vs v
plt.figure(figsize=(12, 8))

plt.subplot(311)
plt.plot(t_points, C_values, label='C')
plt.xlabel('Time')
plt.ylabel('C')
plt.legend()

plt.subplot(312)
plt.plot(t_points, s_values, label='s')
plt.xlabel('Time')
plt.ylabel('s')
plt.legend()

plt.subplot(313)
plt.plot(t_points, v_values, label='v')
plt.xlabel('Time')
plt.ylabel('v')
plt.legend()

plt.tight_layout()
plt.show()

# Print the s values
print("s values:")
for s in s_values:
    print(s)
2

C is just a function of v, so it seems to me you can rewrite the differential equation as

dv / dt = func(v)

and define func in terms of v based on the given conditions?

Perhaps looking at the doc of scipy.integrate.solve_ivp — SciPy v1.11.4 Manual would help? It gives an example where func is an exponential decay function, but in principle you can use any function.

3

Can you provide a screenshot of the equations in written form (non-code form)?
Would probably help to make sure that those viewing the equations are interpreting
them correctly.

dv/dt → acceleration

By the way:

(0.11551-0.002744C v^2)/0.2025 = dv/dt

Can best be written as (to get rid of the constant 0.2025 from the denominator):
dv/dt = 0.57042 - 0.01355Cv^2 (acceleration equation)

Three differential equations for acceleration for different values of C (boundary conditions):

  1. dv/dt = 0.11551 - 0.001322v # Eq. 1 for C = 322/ (3300v) when 3300v<10
  2. Eq. 2 - verify equation stated is correct and plug in for C - can’t really tell if correctly stated - provide screenshot of equation in written form (for 10<3300v<20000)
  3. dv/dt = 0.1551 - 0.010027v^2 # Eq. 3 for C = 0.74 when 3300v>20000

Since you can plug in the C values into the velocity equation for the three different boundary
conditions beforehand, this will give you three differential equations to solve. Since C is a function
of v, this will get rid of C in the differential equations and they will all be a function of v. This
should greatly simplify your code for solving the three differential equations.

1 Like