Programming in python
NGWOH MEH ROY
ICTU20223285
Question 1)
Plot the function
We’ll evaluate f(x) over a range of x values and plot it.
```python
import numpy as np
import [Link] as plt
def f(x):
total = 0.0
for k in range(1, 101):
total += [Link](k**2 * x) / (k**2)
return total
x_vals = [Link](-2*[Link], 2*[Link], 5000)
y_vals = [f(x) for x in x_vals]
[Link](figsize=(10, 5))
[Link](x_vals, y_vals, label=r'$f(x) = \sum_{k=1}^{100} \frac{\sin(k^2 x)}{k^2}$')
[Link]('x')
[Link]('f(x)')
[Link]('Graph of f(x)')
[Link](True)
[Link]()
[Link]()
```
Observation:
The function looks continuous but has many small, rapid oscillations — it’s related to the
Riemann function type, but here k^2 in the argument makes it highly oscillatory for large k.
B) First and second derivatives
Analytically:
f'(x) = \sum_{k=1}^{100} \frac{k^2 \cos(k^2 x)}{k^2} = \sum_{k=1}^{100} \cos(k^2 x)
f''(x) = - \sum_{k=1}^{100} k^2 \sin(k^2 x)
We can compute numerically using finite differences for verification, but let's plot the analytic
derivatives:
```python
def f_prime(x):
return sum([Link](k**2 * x) for k in range(1, 101))
def f_double_prime(x):
return -sum(k**2 * [Link](k**2 * x) for k in range(1, 101))
y_prime = [f_prime(x) for x in x_vals]
y_double_prime = [f_double_prime(x) for x in x_vals]
[Link](figsize=(12, 8))
[Link](3, 1, 1)
[Link](x_vals, y_vals)
[Link]('f(x)')
[Link](True)
[Link](3, 1, 2)
[Link](x_vals, y_prime)
[Link]("f'(x)")
[Link](True)
[Link](3, 1, 3)
[Link](x_vals, y_double_prime)
[Link]("f''(x)")
[Link](True)
plt.tight_layout()
[Link]()
```
Conclusion:
The first derivative f'(x) is the sum of cosine terms with incommensurable frequencies k^2, so it
does not tend to a simple periodic form — it’s highly oscillatory but bounded.
The second derivative f''(x) has amplitude growing with k^2, so it shows wild oscillations — the
function is smooth but has “pathological” derivatives in terms of behavior, possibly nowhere
monotonic or differentiable in the classical sense for infinite sum (Weierstrass-type function).
For finite k=100, it’s still smooth but visually rough.
Exercise 2)
a) Solve \cos x = x
We solve \cos x - x = 0 using [Link] or fsolve.
```python
from [Link] import root
def eq(x):
return [Link](x) - x
sol = root(eq, 0.5) # initial guess 0.5
print(f"Solution to cos(x) = x: x = {sol.x[0]}")
```
Output:
x \approx 0.7390851332151607
B) First formula:
\pi = 4 I, \quad I = \int_0^1 \sqrt{1 - t^2} \, dt
Actually, I is the area of a quarter unit circle, so \pi = 4I.
Let’s compute I numerically.
```python
from [Link] import quad
I, err = quad(lambda t: [Link](1 - t**2), 0, 1)
pi1 = 4 * I
print(f"π estimate from first formula: {pi1}")
```
Second formula:
They wrote:
\pi = 12 \left( I - \frac{\sqrt{3}}{n} \right) \quad\text{with}\quad J = \int_0^{1/2} \sqrt{1 - t^2} \,
dt
But there seems a typo — the second formula is likely from another geometric decomposition.
Let’s check:
Actually J = \int_0^{1/2} \sqrt{1 - t^2} dt is part of a circle segment.
But the given formula \pi = 12(J - \sqrt{3}/n) doesn’t make sense unless n is specified. Possibly
they meant \pi = 12(J - \sqrt{3}/8) from known geometry? Let’s compute J and see if it matches
known identity.
Let’s compute J:
```python
J, errJ = quad(lambda t: [Link](1 - t**2), 0, 0.5)
# Known: J = (pi + sqrt(3)/4 )/6 ? Let's test:
# Actually: ∫√(1-t²)dt from 0 to 1/2 = (sqrt(3)/4 + pi/6)/2? Let's derive properly:
# Let’s skip the confusion — maybe they meant π = 12*(J - sqrt(3)/4) something? But that’s
negative.
# Possibly they meant: π = 12*(J - sqrt(3)/8) ??? Let's check:
# Let’s just compute J and compare:
print(f"J = {J}")
# Known exact: ∫√(1-t²)dt from 0 to 1/2 = (sqrt(3))/8 + pi/12
# So pi = 12*(J - sqrt(3)/8) indeed.
sqrt3_8 = [Link](3)/8
pi2 = 12 * (J - sqrt3_8)
print(f"π estimate from second formula: {pi2}")
```
Both should give \pi \approx 3.14159..
C)
Solve ODEs and plot
1. x'(t) = 1 + t^2 x(t)^2
We use solve_ivp:
```python
from [Link] import solve_ivp
def ode1(t, x):
return 1 + t**2 * x**2
t_span = (0, 1)
x0 = [0]
sol1 = solve_ivp(ode1, t_span, x0, t_eval=[Link](0, 1, 100))
[Link]()
[Link](sol1.t, sol1.y[0])
[Link]("Solution to x'(t) = 1 + t² x(t)²")
[Link]('t')
[Link]('x(t)')
[Link](True)
[Link]()
System:
\frac{dx}{dt} = y, \quad \frac{dy}{dt} = \frac{y}{2} - x - y^3
```python
def ode_system(t, z):
x, y = z
return [y, 0.5*y - x - y**3]
z0 = [1, 0]
t_span = (0, 20)
sol2 = solve_ivp(ode_system, t_span, z0, t_eval=[Link](0, 20, 1000))
[Link]()
[Link](sol2.t, sol2.y[0], label='x(t)')
[Link](sol2.t, sol2.y[1], label='y(t)')
[Link]("Solution to the 2D system")
[Link]('t')
[Link]()
[Link](True)
[Link]()
Bessel-type equation:
t^2 x''(t) + t x'(t) + (t^2 - \alpha^2)x(t) = 0, \quad \alpha=0.5
This is Bessel’s equation; solutions are Bessel functions.
We rewrite as:
x'' = -\frac{1}{t} x' - \left(1 - \frac{\alpha^2}{t^2}\right)x
We can solve numerically with initial conditions e.g. x(0.1) = J_{0.5}(0.1), x'(0.1) = J_{0.5}'(0.1) to
compare.
But for brevity, let’s just solve from t=0.1 to t=10 with approximate ICs from known Bessel
values.
```python
from [Link] import jv, jvp
alpha = 0.5
def bessel_ode(t, z):
x, xp = z
if t == 0:
# avoid singularity, start from small t
xpp = 0
else:
xpp = -xp/t - (1 - alpha**2/t**2)*x
return [xp, xpp]
# Initial conditions at t=0.1
t0 = 0.1
x0 = jv(alpha, t0)
xp0 = jvp(alpha, t0) # derivative of Bessel J
sol3 = solve_ivp(bessel_ode, (t0, 10), [x0, xp0], t_eval=[Link](t0, 10, 300))
[Link]()
[Link](sol3.t, sol3.y[0], label='Numerical')
[Link](sol3.t, jv(alpha, sol3.t), '--', label='Exact Bessel J')
[Link]("Solution to Bessel equation")
[Link]('t')
[Link]('x(t)')
[Link]()
[Link](True)
[Link]()
Exercise 3)