Computer programming can support two related mathematical activities:
A useful Python workflow combines:
The key principle is to choose a computational representation that matches the form of the curve.
An explicit curve is written as
\[ y=f(x). \]
Example:
\[ y=x^3-3x. \]
For each selected value of \(x\), calculate \(y=f(x)\). The resulting points
\[ (x_1,f(x_1)),(x_2,f(x_2)),\ldots,(x_n,f(x_n)) \]
are plotted and joined to create an approximation of the curve.
A parametric curve is described using a parameter \(t\):
\[ x=x(t), \qquad y=y(t). \]
For a circle of radius 1:
\[ x=\cos(t), \qquad y=\sin(t), \qquad 0\leq t\leq 2\pi. \]
Parametric equations are useful for circles, spirals, motion paths, animations, and curves that do not represent \(y\) as a single-valued function of \(x\).
An implicit curve is defined by an equation such as
\[ F(x,y)=0. \]
Example:
\[ x^2+y^2-4=0. \]
This is a circle of radius 2. It is not necessary to solve explicitly for \(y\). Instead, evaluate \(F(x,y)\) over a grid and draw the zero contour.
The following example plots
\[ y=x^3-3x. \]
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 1000)
y = x**3 - 3*x
plt.axhline(0, color="black", linewidth=0.8)
plt.axvline(0, color="black", linewidth=0.8)
plt.plot(x, y, label=r"$y=x^3-3x$")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.legend()
plt.show()
np.linspace(-3, 3, 1000) creates 1,000 equally spaced values between -3 and 3. NumPy then evaluates the expression for all values in the array.
The computer does not draw a mathematical curve directly. It calculates a finite set of points. A larger number of points generally creates a smoother plot, but increasing resolution does not automatically resolve discontinuities, asymptotes, or omitted branches.
Always inspect the domain and special features of the equation before interpreting the plot.
A plot is useful, but calculus explains why the curve has its shape. A standard curve-analysis workflow is:
For
\[ f(x)=x^3-3x, \]
use SymPy as follows:
import sympy as sp
x = sp.symbols("x")
f = x**3 - 3*x
first_derivative = sp.diff(f, x)
second_derivative = sp.diff(f, x, 2)
stationary_points = sp.solve(first_derivative, x)
print(first_derivative)
print(second_derivative)
print(stationary_points)3*x**2 - 3
6*x
[-1, 1]The first derivative is
\[ f'(x)=3x^2-3. \]
The stationary points satisfy
\[ 3x^2-3=0, \]
so \(x=-1\) and \(x=1\). Their corresponding coordinates can be calculated with:
[(p, f.subs(x, p)) for p in stationary_points][(-1, 2), (1, -2)]import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0, 2*np.pi, 1000)
x = np.cos(t)
y = np.sin(t)
plt.plot(x, y)
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.show()
plt.axis("equal") ensures that one unit on the horizontal axis has the same visual length as one unit on the vertical axis. Without it, a circle can appear to be an ellipse.
t = np.linspace(0, 6*np.pi, 2000)
r = 0.1 * t
x = r * np.cos(t)
y = r * np.sin(t)
plt.plot(x, y)
plt.axis("equal")
plt.grid(True)
plt.show()
For
\[ F(x,y)=x^2+y^2-4=0, \]
create a two-dimensional grid and draw the contour where the expression equals zero:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 500)
y = np.linspace(-3, 3, 500)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2 - 4
plt.contour(X, Y, Z, levels=[0])
plt.axis("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.show()
The plotted contour represents all points satisfying
\[ F(x,y)=0. \]
This approach is useful when solving for \(y\) explicitly is difficult or would produce only one branch of the curve.
A root of a function is a value \(x\) for which
\[ f(x)=0. \]
Some equations have exact algebraic solutions. Many important equations do not, so a numerical method is used to obtain an approximation.
brentqConsider
\[ f(x)=x^3-2x-5. \]
import numpy as np
from scipy.optimize import brentq
def f(x):
return x**3 - 2*x - 5
root = brentq(f, 2, 3)
print(root)
print(f(root))2.094551481542327
3.552713678800501e-15The root is approximately
\[ x\approx 2.0945514815. \]
The interval \([2,3]\) is suitable because the function changes sign across the interval:
f(2), f(3)(-1, 16)A sign change indicates that a root is bracketed, assuming the function is continuous on the interval. brentq is a reliable bracketing method for one-dimensional roots.
import matplotlib.pyplot as plt
x = np.linspace(-3, 3, 500)
y = f(x)
plt.axhline(0, color="black", linewidth=0.8)
plt.plot(x, y, label=r"$f(x)=x^3-2x-5$")
plt.scatter(root, 0, color="red", zorder=3,
label=f"root ≈ {root:.3f}")
plt.grid(True)
plt.legend()
plt.show()
If the intervals containing roots are unknown, sample the function over a range and look for sign changes.
from scipy.optimize import brentq
import numpy as np
def find_roots_by_scan(function, left, right, points=10000):
xs = np.linspace(left, right, points)
ys = function(xs)
roots = []
for x1, x2, y1, y2 in zip(xs[:-1], xs[1:], ys[:-1], ys[1:]):
if y1 == 0:
roots.append(x1)
elif y1 * y2 < 0:
roots.append(brentq(function, x1, x2))
return np.unique(np.round(roots, 12))Example:
def f(x):
return np.sin(x)
roots = find_roots_by_scan(f, -2*np.pi, 2*np.pi)
print(roots)[-3.14159265 0. 3.14159265]A sign-change scan can miss a root where the curve touches the horizontal axis without crossing it. For example:
\[ f(x)=(x-2)^2. \]
The function has a root at \(x=2\), but it is positive on both sides of the root. Consequently, there is no sign change.
Possible responses include:
nsolve.Suppose two curves are
\[ y=f(x), \qquad y=g(x). \]
At an intersection,
\[ f(x)=g(x). \]
Rearrange this as
\[ h(x)=f(x)-g(x)=0. \]
Thus, curve intersections are roots of the difference function.
Let
\[ f(x)=x^2, \qquad g(x)=2x+3. \]
def f(x):
return x**2
def g(x):
return 2*x + 3
def difference(x):
return f(x) - g(x)Find the two roots of the difference:
x1 = brentq(difference, -2, 0)
x2 = brentq(difference, 3, 5)
intersections = [
(x1, f(x1)),
(x2, f(x2))
]
print(intersections)[(-0.9999999999999986, 0.9999999999999971), (3.0, 9.0)]The intersections are
\[ (-1,1) \quad\text{and}\quad (3,9). \]
Plot the curves and mark the points:
x = np.linspace(-2, 5, 500)
plt.plot(x, f(x), label=r"$y=x^2$")
plt.plot(x, g(x), label=r"$y=2x+3$")
for x_value, y_value in intersections:
plt.scatter(x_value, y_value, color="red", zorder=3)
plt.axhline(0, color="black", linewidth=0.8)
plt.axvline(0, color="black", linewidth=0.8)
plt.grid(True)
plt.legend()
plt.show()
def find_intersections(f, g, left, right, points=5000):
h = lambda x: f(x) - g(x)
xs = np.linspace(left, right, points)
ys = h(xs)
roots = []
for a, b, fa, fb in zip(xs[:-1], xs[1:], ys[:-1], ys[1:]):
if fa * fb < 0:
roots.append(brentq(h, a, b))
roots = np.unique(np.round(roots, 10))
return [(x_value, f(x_value)) for x_value in roots]Use the function as follows:
f = lambda x: x**2
g = lambda x: 2*x + 3
points = find_intersections(f, g, -5, 5)
print(points)[(np.float64(-1.0), np.float64(1.0)), (np.float64(3.0), np.float64(9.0))]The algorithm is:
brentq to each bracket.SymPy is helpful when the equations are represented symbolically.
import sympy as sp
x = sp.symbols("x")
equation = x**3 - 2*x - 5
root = sp.nsolve(equation, 2)
print(root)2.09455148154233The second argument is an initial guess. For the intersection example:
f_expr = x**2
g_expr = 2*x + 3
intersection_x = sp.nsolve(f_expr - g_expr, -1)
intersection_y = f_expr.subs(x, intersection_x)
print(intersection_x)
print(intersection_y)-1.00000000000000
1.00000000000000Use another initial guess to find the other intersection:
intersection_x = sp.nsolve(f_expr - g_expr, 3)
intersection_y = f_expr.subs(x, intersection_x)
print(intersection_x)
print(intersection_y)3.00000000000000
9.00000000000000nsolve may converge to different roots depending on the initial guess. It is therefore good practice to combine it with a plot, an algebraic analysis, or a systematic scan.
For two implicit curves, solve two equations in two unknowns. Consider:
\[ x^2+y^2=25, \]
and
\[ y=x+1. \]
x, y = sp.symbols("x y")
equations = [
x**2 + y**2 - 25,
y - x - 1
]
solution = sp.nsolve(equations, [x, y], [3, 4])
print(solution)Matrix([[3.00000000000000], [4.00000000000000]])The values [3, 4] are an initial guess for the point \((x,y)\). Try a different initial guess to search for another solution:
solution = sp.nsolve(equations, [x, y], [-4, -3])
print(solution)Matrix([[-4.00000000000000], [-3.00000000000000]])For multidimensional numerical root-finding, SciPy also provides scipy.optimize.root.
Newton’s method uses the recurrence
\[ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}. \]
A simple implementation is:
def newton(function, derivative, initial_guess,
tolerance=1e-10, max_iterations=100):
x = initial_guess
for _ in range(max_iterations):
next_x = x - function(x) / derivative(x)
if abs(next_x - x) < tolerance:
return next_x
x = next_x
raise RuntimeError("Newton's method did not converge")Example:
def f(x):
return x**3 - 2*x - 5
def derivative(x):
return 3*x**2 - 2
root = newton(f, derivative, initial_guess=2)
print(root)2.0945514815423265Newton’s method is fast near a root, but it can fail when the initial guess is poor, the derivative is zero or very small, or the iteration moves into an unsuitable region. For robust general-purpose work, a bracketing method such as Brent’s method is often preferable when a valid interval is available.
| Task | Recommended approach |
|---|---|
| Plot \(y=f(x)\) | NumPy and Matplotlib |
| Plot a circle, spiral, or motion path | Parametric equations |
| Plot \(F(x,y)=0\) | Matplotlib contour |
| Differentiate or integrate | SymPy |
| Find exact algebraic roots | SymPy solve |
| Find one bracketed numerical root | SciPy brentq |
| Find roots when intervals are unknown | Scan and then use brentq |
| Find a root from an initial guess | SymPy nsolve or Newton’s method |
| Solve several equations in several variables | SymPy nsolve or SciPy root |
Before accepting a result:
brentq, verify that the endpoints have opposite signs.abs(f(root)).nsolve or Newton’s method.The following program finds and plots the intersections between
\[ y=\sin(x) \]
and
\[ y=0.3x. \]
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
f = lambda x: np.sin(x)
g = lambda x: 0.3 * x
h = lambda x: f(x) - g(x)
xs = np.linspace(-10, 10, 10000)
ys = h(xs)
roots = []
for a, b, fa, fb in zip(xs[:-1], xs[1:], ys[:-1], ys[1:]):
if fa == 0:
roots.append(a)
elif fa * fb < 0:
roots.append(brentq(h, a, b))
roots = np.unique(np.round(roots, 10))
intersection_points = [(r, f(r)) for r in roots]
plt.plot(xs, f(xs), label=r"$y=\sin(x)$")
plt.plot(xs, g(xs), label=r"$y=0.3x$")
for x_value, y_value in intersection_points:
plt.scatter(x_value, y_value, color="red", zorder=3)
plt.axhline(0, color="black", linewidth=0.8)
plt.axvline(0, color="black", linewidth=0.8)
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.legend()
plt.show()
print(intersection_points)
[(np.float64(-2.3564411499), np.float64(-0.7069323449258418)), (np.float64(0.0), np.float64(0.0)), (np.float64(2.3564411499), np.float64(0.7069323449258418))]This example illustrates the complete pattern:
\[ \text{two curves} \rightarrow \text{difference function} \rightarrow \text{roots} \rightarrow \text{intersection coordinates} \rightarrow \text{visual verification}. \]
nsolve with different initial guesses and record which solution is found.brentq.The most reusable idea in computational curve solving is:
\[ \boxed{\text{An intersection of } f \text{ and } g \text{ is a root of } f-g.} \]
Use plots to explore the geometry, symbolic tools to analyse the mathematics, and numerical solvers to obtain approximate values when exact solutions are unavailable.