Advanced Calculus with SageMath - Chapter 9

# Comprehensive imports for advanced calculus
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.optimize as opt
import scipy.integrate as integrate
from scipy.integrate import solve_ivp, odeint
import sympy as sp
from sympy import *
from sage.all import *
import seaborn as sns

# Configure plotting
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")
plt.rcParams['figure.figsize'] = (12, 8)

print("Advanced Calculus Environment Initialized")
print("Tools: SageMath, NumPy, SciPy, SymPy, Matplotlib")
print("Ready for multivariable calculus, vector analysis, and PDEs!")

# Advanced computational techniques
print("ADVANCED COMPUTATIONAL METHODS")
print("=" * 50)

# Automatic differentiation example using optimization
from scipy.optimize import minimize

# Define a complex multivariable function
def complex_function(vars):
    x, y = vars
    return np.sin(x*y) * np.exp(-(x**2 + y**2)) + 0.1*(x**4 + y**4)

print("Complex function: f(x,y) = sin(xy)·exp(-(x²+y²)) + 0.1(x⁴+y⁴)")

# Find minimum using optimization
result = minimize(complex_function, [1.0, 1.0], method='BFGS')

print("\nOPTIMIZATION RESULT")
print(f"Minimum found at: x = {result.x[0]:.6f}, y = {result.x[1]:.6f}")
print(f"Function value: f = {result.fun:.6f}")
print(f"Convergence: {result.success}")
print(f"Iterations: {result.nit}")

# Monte Carlo integration example
print("\nMONTE CARLO INTEGRATION")
print("Estimating ∬_D e^(-(x²+y²)) dx dy over unit disk")

def monte_carlo_integral(n_samples=100000):
    # Generate random points in unit square
    x_rand = np.random.uniform(-1, 1, n_samples)
    y_rand = np.random.uniform(-1, 1, n_samples)
    
    # Check which points are inside unit disk
    inside_disk = (x_rand**2 + y_rand**2) <= 1
    
    # Evaluate function at points inside disk
    f_values = np.exp(-(x_rand[inside_disk]**2 + y_rand[inside_disk]**2))
    
    # Estimate integral
    disk_area = np.pi
    square_area = 4
    n_inside = np.sum(inside_disk)
    
    integral_estimate = (square_area * np.sum(f_values) / n_inside) * (n_inside / n_samples)
    
    return integral_estimate, n_inside

mc_result, points_inside = monte_carlo_integral()
analytical_result = np.pi * (1 - np.exp(-1))  # Analytical answer

print(f"Monte Carlo estimate: {mc_result:.6f}")
print(f"Analytical result: {analytical_result:.6f}")
print(f"Error: {abs(mc_result - analytical_result):.6f}")
print(f"Points inside disk: {points_inside}")

# Symbolic computation with coordinate transformations
print("\nSYMBOLIC COMPUTATION")
print("Vector field transformations and calculus operations")

# Vector field in cylindrical coordinates
rho, phi = symbols('rho phi', positive=True)
F_rho = rho * cos(phi)
F_phi = sin(phi)
F_z_cyl = rho**2

print(f"Vector field in cylindrical coordinates:")
print(f"F_ρ = {F_rho}")
print(f"F_φ = {F_phi}")
print(f"F_z = {F_z_cyl}")

# This demonstrates the power of symbolic computation
print("\nThis demonstrates the power of symbolic computation")
print("for complex coordinate transformations and calculations!")

import numpy as np
import  as plt

# Illustrations: optimization landscape, Monte Carlo integration, and a vector field projection

# 1) Contour of the complex function with the BFGS minimum
x = np.linspace(-2, 2, 300)
y = np.linspace(-2, 2, 300)
X, Y = np.meshgrid(x, y)
Z = np.sin(X*Y) * np.exp(-(X**2 + Y**2)) + 0.1*(X**4 + Y**4)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

cs = axes[0].contourf(X, Y, Z, levels=40, cmap='viridis')
fig.colorbar(cs, ax=axes[0], shrink=0.9)
axes[0].set_title('Complex function f(x, y)')
axes[0].set_xlabel('x'); axes[0].set_ylabel('y')
try:
    xm, ym = result.x
    axes[0].plot(xm, ym, 'r*', ms=12, label='BFGS minimum')
    axes[0].legend(loc='upper right')
except Exception:
    pass

# 2) Monte Carlo sampling in the unit disk for ∫∫ e^{-(x^2+y^2)} dx dy
np.random.seed(0)
N = 5000
xr = np.random.uniform(-1, 1, N)
yr = np.random.uniform(-1, 1, N)
inside = (xr**2 + yr**2) <= 1.0
fvals = np.exp(-(xr**2 + yr**2))
mc_est = 4.0 * np.sum(fvals * inside) / N
analytical = np.pi * (1 - np.exp(-1))

axes[1].scatter(xr[~inside], yr[~inside], s=3, color='0.85', alpha=0.3, label='Outside')
p = axes[1].scatter(xr[inside], yr[inside], s=6, c=fvals[inside], cmap='plasma', alpha=0.8, label='Inside')
theta = np.linspace(0, 2*np.pi, 400)
axes[1].plot(np.cos(theta), np.sin(theta), 'k--', lw=1)
fig.colorbar(p, ax=axes[1], label='e^{-(x^2+y^2)}', shrink=0.9)
axes[1].set_aspect('equal', 'box')
axes[1].set_xlim(-1.1, 1.1); axes[1].set_ylim(-1.1, 1.1)
axes[1].set_title('Monte Carlo in the unit disk')
axes[1].text(-1.05, -1.3, f"MC ≈ {mc_est:.4f}\nAnalytical ≈ {analytical:.4f}", fontsize=9)
axes[1].legend(loc='upper right')

# 3) Vector field from cylindrical coordinates projected to the xy-plane
gx = np.linspace(-2, 2, 21)
gy = np.linspace(-2, 2, 21)
GX, GY = np.meshgrid(gx, gy)
R = np.sqrt(GX**2 + GY**2)
Phi = np.arctan2(GY, GX)

Fr = R * np.cos(Phi)
Fphi = np.sin(Phi)

Fx = Fr * np.cos(Phi) - Fphi * np.sin(Phi)
Fy = Fr * np.sin(Phi) + Fphi * np.cos(Phi)

Fx = np.nan_to_num(Fx)
Fy = np.nan_to_num(Fy)
mag = np.sqrt(Fx**2 + Fy**2); mag[mag == 0] = 1.0
Fxn = Fx / mag
Fyn = Fy / mag

axes[2].quiver(GX, GY, Fxn, Fyn, color='teal', pivot='mid', angles='xy', scale_units='xy', scale=0.9)
axes[2].set_title('Vector field (cylindrical → xy projection)')
axes[2].set_xlim(-2, 2); axes[2].set_ylim(-2, 2)
axes[2].set_aspect('equal', 'box')
axes[2].set_xlabel('x'); axes[2].set_ylabel('y')

plt.tight_layout()
plt.show()