Exercises — SymPy — symbolic algebra, calculus, ODE solving

WHAT: Opening brackets = multiplying everything out. WHY that tool: expand is the "multiply out" pencil move; factor is its reverse, simplify is a heuristic that might or might not do it.

x = symbols('x')
expand((x + 3)*(x - 3))     # -> x**2 - 9

Result: x**2 - 9.

WHAT: finding roots. WHY solve: solve(eq, x) treats eq as eq = 0 and returns a list of symbolic roots — not a single number.

x = symbols('x')
solve(x**2 - 7*x + 12, x)    # -> [3, 4]

Result: [3, 4]. Note it is a list, so you would index sol[0] to grab 3.

WHAT: swap the symbol, keep exactness. WHY .subs then simplify: .subs swaps x for sqrt(3) but leaves the expression un-collapsed; simplify combines it. What it looks like: $(\sqrt{3}-2{)}^2=3-4\sqrt{3}+4=7-4\sqrt{3}$.

x = symbols('x')
e = (x - 2)**2
simplify(e.subs(x, sqrt(3)))   # -> 7 - 4*sqrt(3)

Result: 7 - 4*sqrt(3). If you wanted a decimal you'd add .evalf() → 0.0718....

WHAT: differentiate a product. WHY diff: it applies the known rules (product + power + trig) automatically.

x = symbols('x')
diff(x**3 * sin(x), x)     # -> x**3*cos(x) + 3*x**2*sin(x)

Result: x**3*cos(x) + 3*x**2*sin(x). The two rules: the product rule $(uv{)}^{\prime }=u^{\prime }v+u{v}^{\prime }$ with $u=x^3$ (power rule gives $u^{\prime }=3x^2$) and $v=\sin x$ (gives $v^{\prime }=\cos x$).

WHAT: area under the curve on $[0,2]$. WHY the tuple form (x,a,b): it applies the Fundamental Theorem ${\int }_a^{b}f=F(b)-F(a)$. By hand: $F=x^3+x$, so $F(2)-F(0)=(8+2)-0=10$.

x = symbols('x')
integrate(3*x**2 + 1, (x, 0, 2))   # -> 10

Result: 10 (an exact integer, not 10.0).

WHAT: expand $\cos x$ around $0$ up to (but not including) order 7.

x = symbols('x')
cos(x).series(x, 0, 7)
# -> 1 - x**2/2 + x**4/24 - x**6/720 + O(x**7)

Result: 1 - x**2/2 + x**4/24 - x**6/720 + O(x**7). WHY only even powers: the Taylor coefficient of $x^n$ is $\frac{f^{(n)}(0)}{n!}$ (see Taylor and Maclaurin Series). For $\cos$, every odd derivative at $0$ is a $\pm \sin (0)=0$, so those coefficients vanish. Geometrically, $\cos$ is an even function (mirror-symmetric across the $y$-axis), and even functions have only even-power terms. Look at figure s01: the mirror symmetry is the visual reason.

WHAT: roots of a quadratic with no real solutions.

x = symbols('x')
solve(x**2 + x + 1, x)
# -> [-1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]

Result: [-1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]. WHY this form: the quadratic formula gives $x=\frac{-1\pm \sqrt{1-4}}{2}=\frac{-1\pm \sqrt{-3}}{2}$. Since $\sqrt{-3}=\sqrt{3}i$ (where $i^2=-1$), SymPy keeps it exact and complex rather than crashing. This is the same complex-root logic the characteristic equation method uses to turn $r=\pm i$ into oscillations.

WHAT: compute ${\lim }_{h\to 0}\frac{(x+h{)}^3-x^3}{h}$. WHY the limit and not diff: we want to derive, not assume — this is the definition from Limits and the definition of the derivative.

x, h = symbols('x h')
limit(((x + h)**3 - x**3)/h, h, 0)   # -> 3*x**2
diff(x**3, x)                        # -> 3*x**2  (matches!)

Result: both give 3*x**2. What it looks like: expanding, $(x+h{)}^3-x^3=3x^{2}h+3x{h}^2+h^3$; divide by $h$ to get $3x^2+3xh+h^2$; let $h\to 0$ and only $3x^2$ survives.

WHAT: three tools chained — series, .removeO(), and lambdify. WHY chain them: series gives a symbolic polynomial; .removeO() strips the bookkeeping $O(x^6)$ so it becomes a real polynomial; lambdify compiles it into a fast numeric function for evaluation.

x = symbols('x')
p = sin(x).series(x, 0, 6).removeO()   # x - x**3/6 + x**5/120
f = lambdify(x, p, 'math')
f(0.5)          # -> 0.4794270833333333
math.sin(0.5)   # -> 0.479425538604203

The polynomial gives 0.47942708..., the true value is 0.47942553... — they agree to 5 decimals because we kept terms through $x^5$. What it looks like: figure s02 overlays the polynomial (orange) on true $\sin$ (blue); near 0 they hug, far away they peel apart.

FORECAST: $y^{\prime }=3y$ means "growth rate proportional to size" → exponential $e^{3x}$. General solution $C_1{e}^{3x}$; $y(0)=5$ pins $C_1=5$. So predict $y=5e^{3x}$.

x = symbols('x')
y = Function('y')
sol = dsolve(Eq(y(x).diff(x) - 3*y(x), 0), y(x), ics={y(0): 5})
sol            # -> Eq(y(x), 5*exp(3*x))

Result: Eq(y(x), 5*exp(3*x)) — matches the forecast. Verify: $y^{\prime }=15e^{3x}$ and $3y=15e^{3x}$, so $y^{\prime }-3y=0$ ✓ and $y(0)=5e^0=5$ ✓.

WHAT: distinguish structural identity from mathematical equality.

x = symbols('x')
(x + 2)**2 == x**2 + 4*x + 4        # -> False  (different trees!)
simplify((x + 2)**2 - (x**2 + 4*x + 4)) == 0   # -> True
((x + 2)**2).equals(x**2 + 4*x + 4)            # -> True

WHY False: == compares the shape of the expression tree. (x+2)**2 is stored as a power node; x**2+4*x+4 is a sum node — different structures, so == is False even though they're equal as functions. Fix: subtract and simplify to 0, or use .equals.

WHAT: improper integrals (upper limit is oo, SymPy's infinity).

x = symbols('x')
integrate(1/x**2, (x, 1, oo))   # -> 1
integrate(1/x,   (x, 1, oo))    # -> oo

Results: 1 and oo. WHY the split: the antiderivative of $1/x^2$ is $-1/x$, which approaches $0$ as $x\to \infty$, so the area is bounded: $[-1/x{]}_1^{\infty }=0-(-1)=1$. But the antiderivative of $1/x$ is $\ln x$, which grows without bound, so the area is infinite. Edge lesson: two curves that both shrink toward zero can enclose wildly different total areas — $1/x^2$ shrinks fast enough, $1/x$ does not.

WHAT: simplification depends on what SymPy is allowed to assume.

x = symbols('x')
simplify(sqrt(x**2))          # -> sqrt(x**2)  (stays, or Abs(x))
p = symbols('p', positive=True)
simplify(sqrt(p**2))          # -> p

WHY: for a general real $x$, $\sqrt{x^2}=|x|$ (think $x=-3$: $\sqrt{9}=3=|-3|$, not $-3$). SymPy refuses to drop the absolute value because it can't prove $x\ge 0$. Once you declare positive=True, the assumption lets it simplify to just p. Edge lesson: the same expression has two correct answers depending on the sign domain — always tell SymPy your assumptions.

WHAT: combine solve, list indexing, and .evalf.

x = symbols('x')
roots = solve(x**2 - 2, x)     # -> [-sqrt(2), sqrt(2)]
pos = roots[1]                 # sqrt(2)
pos                            # exact: sqrt(2)
pos.evalf(6)                   # -> 1.41421

Results: exact root sqrt(2), decimal 1.41421. WHY index then evalf: solve returns a list of exact symbolic roots; we pick the positive one, keep it exact for algebra, and only collapse to a float at the very end.

(x+2)**2 == x**2+4*x+4 return False? == compares expression-tree structure, not mathematical value. Use simplify(a-b)==0 or a.equals(b).

lambdify-ing a series? .removeO(), to strip the O(x**n) order term that isn't numerically evaluable.

${\int }_1^{\infty }\frac{1}{x}dx$ infinite but ${\int }_1^{\infty }\frac{1}{x^2}dx$ finite? The antiderivative $\ln x$ grows without bound; $-1/x$ approaches a finite limit ($1$).

sqrt(x**2) to x? Because for general real $x$ it equals $|x|$. Declare positive=True to allow the reduction.