a3q1
In [1]: import numpy as np
import [Link] as plt
from [Link] import odeint, solve_ivp
# -------------------------------
# Problem (i)
# -------------------------------
def f1(y, t):
return t*[Link](3*t) - 2*y
def f1_ivp(t, y):
return t*[Link](3*t) - 2*y
t1 = [Link](0, 1, 100)
y0_1 = 0
# odeint
y_odeint_1 = odeint(f1, y0_1, t1)
# solve_ivp
sol1 = solve_ivp(f1_ivp, [0,1], [y0_1], t_eval=t1)
y_ivp_1 = sol1.y[0]
# exact solution
y_exact_1 = (1/5)*t1*[Link](3*t1) - (1/25)*[Link](3*t1) + (1/25)*[Link](-2*t1)
# error
error_odeint_1 = [Link](y_exact_1 - y_odeint_1.flatten())
error_ivp_1 = [Link](y_exact_1 - y_ivp_1)
[Link]()
[Link](t1, y_exact_1, 'k', label="Exact")
[Link](t1, y_odeint_1, '--', label="odeint")
[Link](t1, y_ivp_1, ':', label="solve_ivp")
[Link]()
[Link]("Problem 1(i)")
[Link]()
# -------------------------------
# Problem (ii)
# -------------------------------
def f2(y, t):
return 1 + (t - y)**2
def f2_ivp(t, y):
return 1 + (t - y)**2
t2 = [Link](2, 3, 100)
y0_2 = 1
y_odeint_2 = odeint(f2, y0_2, t2)
sol2 = solve_ivp(f2_ivp, [2,3], [y0_2], t_eval=t2)
y_ivp_2 = sol2.y[0]
# exact solution
y_exact_2 = t2 + 1/(1 - t2)
error_odeint_2 = [Link](y_exact_2 - y_odeint_2.flatten())
error_ivp_2 = [Link](y_exact_2 - y_ivp_2)
[Link]()
[Link](t2, y_exact_2, 'k', label="Exact")
[Link](t2, y_odeint_2, '--', label="odeint")
[Link](t2, y_ivp_2, ':', label="solve_ivp")
[Link]()
[Link]("Problem 1(ii)")
[Link]()
a3q2
In [2]: import numpy as np
import [Link] as plt
from [Link] import solve_ivp
# Define the system
def lotka_volterra(t, z):
x, y = z
dxdt = -0.1*x + 0.02*x*y
dydt = 0.2*y - 0.025*x*y
return [dxdt, dydt]
# Initial conditions
x0 = 6
y0 = 6
z0 = [x0, y0]
# Time span
t_span = (0, 200)
t_eval = [Link](0, 200, 1000)
# Solve using solve_ivp
sol = solve_ivp(lotka_volterra, t_span, z0, t_eval=t_eval)
# Extract solutions
t = sol.t
x = sol.y[0]
y = sol.y[1]
# Plot populations vs time
[Link]()
[Link](t, x, label='Predators (x)')
[Link](t, y, label='Prey (y)')
[Link]('Time')
[Link]('Population (thousands)')
[Link]('Lotka-Volterra Model')
[Link]()
[Link]()
[Link]()
# Find first time when x ≈ y
diff = [Link](x - y)
idx = [Link](diff < 0.01)[0]
if len(idx) > 1:
t_equal = t[idx[1]] # skip t=0
print(f"Approximate time when populations first equal: t ≈ {t_equal:.3f}")
else:
print("No intersection found in given time range.")
Approximate time when populations first equal: t ≈ 135.335
a2q3
In [3]: import numpy as np
import [Link] as plt
from [Link] import solve_ivp
# -------------------------------
# Competition model
# -------------------------------
def model(t, z):
x, y = z
dxdt = x*(2 - 0.4*x - 0.3*y)
dydt = y*(1 - 0.1*y - 0.3*x)
return [dxdt, dydt]
cases = [
[1.5, 3.5],
[1, 1],
[2, 7],
[4.5, 0.5]
]
t_eval = [Link](0, 50, 1000)
# -------------------------------
# Create ONE figure with 4 subplots
# -------------------------------
fig, axes = [Link](2, 2, figsize=(12, 8))
for i, z0 in enumerate(cases):
sol = solve_ivp(model, [0, 50], z0, t_eval=t_eval)
t = sol.t
x = sol.y[0]
y = sol.y[1]
ax = axes[i // 2, i % 2]
[Link](t, x, 'b', label='x(t)')
[Link](t, y, 'r', label='y(t)')
ax.set_title(f'Case {i+1}: x0={z0[0]}, y0={z0[1]}')
ax.set_xlabel('Time')
ax.set_ylabel('Population')
[Link](True, alpha=0.3)
[Link]()
[Link]('Competition Model – All Cases', fontsize=14, fontweight='bold')
plt.tight_layout()
[Link]()
a3q4
In [4]: import numpy as np
import pandas as pd
import [Link] as plt
from [Link] import solve_ivp
g = 32.17 # ft/s²
L = 2.0 # ft
def pendulum(t, z):
theta, omega = z
dtheta_dt = omega
domega_dt = -(g / L) * [Link](theta)
return [dtheta_dt, domega_dt]
theta0 = [Link] / 6 # π/6 radians ≈ 30°
omega0 = 0.0
t_span = (0, 2)
h = 0.1
t_eval = [Link](0, 2 + h, h) # 0, 0.1, 0.2, ..., 2.0
sol = solve_ivp(pendulum, t_span, [theta0, omega0], t_eval=t_eval)
t = sol.t
theta = sol.y[0]
omega = sol.y[1]
data = [Link]({
"Time (s)": t,
"Theta (rad)": theta,
"Omega (rad/s)": omega
})
# Show the full table (or you can use head()/tail())
print(data)
# -------------------------
# Plot
# -------------------------
[Link]()
[Link](t, theta, label="Theta (rad)")
[Link](t, omega, label="Omega (rad/s)")
[Link]("Time (s)")
[Link]("Values")
[Link]("Pendulum Motion")
[Link]()
[Link]()
[Link]()
Time (s) Theta (rad) Omega (rad/s)
0 0.0 0.523599 0.000000
1 0.1 0.483854 -0.785591
2 0.2 0.370229 -1.460132
3 0.3 0.199187 -1.917018
4 0.4 -0.003293 -2.076131
5 0.5 -0.205173 -1.907780
6 0.6 -0.374796 -1.440051
7 0.7 -0.486320 -0.760615
8 0.8 -0.523396 0.027172
9 0.9 -0.480898 0.811263
10 1.0 -0.365067 1.479569
11 1.1 -0.192304 1.926839
12 1.2 0.010535 2.075615
13 1.3 0.211690 1.894580
14 1.4 0.379650 1.419010
15 1.5 0.488745 0.734399
16 1.6 0.522965 -0.055397
17 1.7 0.478040 -0.836798
18 1.8 0.359712 -1.498472
19 1.9 0.185532 -1.936713
20 2.0 -0.017676 -2.073106
a3q5
In [5]: import numpy as np
from [Link] import solve_ivp
import [Link] as plt
def system(t, u):
u1, u2, u3, u4, u5, u6 = u
du1 = u2
du2 = u3
du3 = -2.0 * u5**2 + u4 # x1''' = -2*(x2')^2 + x2
du4 = u5
du5 = u6
du6 = -u3**3 + u5 + u1 + [Link](t) # x2''' = -x1'' + x2' + x1 + sin(t)
return [du1, du2, du3, du4, du5, du6]
t_span = (0, 10)
t_eval = [Link](0, 10, 500)
u0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
solution = solve_ivp(system, t_span, u0, t_eval=t_eval, method='RK45')
# 4. Extract and plot the results
t = solution.t
x1 = solution.y[0] # u1 corresponds to x1
x2 = solution.y[3] # u4 corresponds to x2
[Link](figsize=(10, 6))
[Link](t, x1, label='$x_1(t)$', color='blue')
[Link](t, x2, label='$x_2(t)$', color='red')
[Link]("Numerical Solution of the System of ODEs")
[Link]("Time (t)")
[Link]("State Variables")
[Link](True)
[Link]()
plt.tight_layout()
[Link]()
a3q6
In [6]: import numpy as np
import pandas as pd
import [Link] as plt
from [Link] import solve_ivp, solve_bvp
# ── Exact solution ────────────────────────────────────────────
def exact(x):
return [Link](-10*x)
# ── Linear Shooting Algorithm ─────────────────────────────────
def ode(x, z):
return [z[1], 100*z[0]]
# Linear shooting (fixed accuracy)
def shooting(h):
x = [Link](0, 1+h, h)
ya, yb = 1, [Link](-10)
u = solve_ivp(ode, [0,1], [ya,0], t_eval=x).y[0]
v = solve_ivp(ode, [0,1], [0,1], t_eval=x).y[0]
t = (yb - u[-1]) / v[-1]
y = u + t*v
return x, y
# solve_bvp (fixed initial guess)
def fun(x, y):
return [Link]([y[1], 100*y[0]])
def bc(ya, yb):
return [Link]([ya[0]-1, yb[0]-[Link](-10)])
x_init = [Link](0,1,10)
y_init = [Link]((2, x_init.size))
y_init[0] = exact(x_init) # IMPORTANT FIX
sol = solve_bvp(fun, bc, x_init, y_init, tol=1e-8)
# Tables
for h in [0.1, 0.05]:
x, y = shooting(h)
df = [Link]({
"x": x,
"Shooting": y,
"Exact": exact(x),
"Error": [Link](y - exact(x))
})
print(f"\nTable (h = {h})")
print(df.to_string(index=False))
# Plot
x_fine = [Link](0,1,400)
[Link](x_fine, exact(x_fine), label="Exact")
x1, y1 = shooting(0.1)
x2, y2 = shooting(0.05)
[Link](x1, y1, 'o--', label="h=0.1")
[Link](x2, y2, 's--', label="h=0.05")
x_bvp = [Link](0,1,100)
[Link](x_bvp, [Link](x_bvp)[0], label="solve_bvp")
[Link]("x")
[Link]("y")
[Link]("BVP Solution")
[Link]()
[Link]()
[Link]()
Table (h = 0.1)
x Shooting Exact Error
0.0 1.000000 1.000000 0.000000e+00
0.1 0.367698 0.367879 1.814879e-04
0.2 0.135229 0.135335 1.058205e-04
0.3 0.050206 0.049787 4.194110e-04
0.4 0.019143 0.018316 8.272164e-04
0.5 0.004969 0.006738 1.769436e-03
0.6 -0.001546 0.002479 4.024734e-03
0.7 0.008517 0.000912 7.604680e-03
0.8 0.016145 0.000335 1.581004e-02
0.9 -0.179803 0.000123 1.799260e-01
1.0 0.000045 0.000045 6.130593e-13
Table (h = 0.05)
x Shooting Exact Error
0.00 1.000000 1.000000 0.000000e+00
0.05 0.606483 0.606531 4.753113e-05
0.10 0.367698 0.367879 1.814879e-04
0.15 0.223127 0.223130 2.817216e-06
0.20 0.135229 0.135335 1.058205e-04
0.25 0.081949 0.082085 1.356428e-04
0.30 0.050206 0.049787 4.194110e-04
0.35 0.029669 0.030197 5.278946e-04
0.40 0.019143 0.018316 8.272164e-04
0.45 0.011499 0.011109 3.898124e-04
0.50 0.004969 0.006738 1.769436e-03
0.55 0.007699 0.004087 3.612066e-03
0.60 -0.001546 0.002479 4.024734e-03
0.65 -0.000946 0.001503 2.449354e-03
0.70 0.008517 0.000912 7.604680e-03
0.75 -0.033449 0.000553 3.400233e-02
0.80 0.016145 0.000335 1.581004e-02
0.85 0.042056 0.000203 4.185270e-02
0.90 -0.179803 0.000123 1.799260e-01
0.95 0.169547 0.000075 1.694717e-01
1.00 0.000045 0.000045 6.130593e-13
a3q7
In [7]: import numpy as np
from [Link] import solve_ivp
import [Link] as plt
def rlc_system(t, y, R, L, C):
q, i = y
V_t = 0.0 # No external voltage
dq_dt = i
di_dt = (1/L) * (V_t - R*i - q/C)
return [dq_dt, di_dt]
# 2. Setup Parameters and Initial Conditions
L = 1.0
C = 1.0
t_span = (0, 20)
t_eval = [Link](0, 20, 1000)
# Initial conditions: q(0) = 1, i(0) = 0
y0 = [1.0, 0.0]
# List of R values to test (starting with 0.5, then reducing)
R_values = [0.5, 0.1, 0.0]
colors = ['blue', 'orange', 'green']
labels = ['R = 0.5 (Original)', 'R = 0.1 (Reduced)', 'R = 0.0 (No resistance)']
# 3. Solve and Plot
[Link](figsize=(10, 6))
for R, color, label in zip(R_values, colors, labels):
# Solve the IVP for the current R value
solution = solve_ivp(rlc_system, t_span, y0, args=(R, L, C), t_eval=t_eval)
# Extract charge q(t)
t = solution.t
q = solution.y[0]
# Plot q(t)
[Link](t, q, color=color, label=label, linewidth=2)
# Styling the plot
[Link]('RLC Circuit: Charge $q(t)$ over Time', fontsize=14)
[Link]('Time ($t$)', fontsize=12)
[Link]('Charge ($q$)', fontsize=12)
[Link](0, color='black', linewidth=0.7) # Add x-axis line
[Link](True, linestyle='--', alpha=0.7)
[Link](loc='upper right')
[Link](0, 20)
plt.tight_layout()
[Link]()