# Import plotting functions
import matplotlib.pyplot as plt

# Sequence an = 8n!
an = [8 * factorial(n) for n in range(1, 11)]
an_conclusion = "The sequence an represents 8n! and converges to 0 as n approaches infinity. It is a monotone decreasing sequence."

# Sequence bn = n + 1
bn = [n + 1 for n in range(1, 11)]
bn_conclusion = "The sequence bn is defined as n + 1 and does not converge. It is not a monotone sequence."

# Sequence cn = [ln(n)]^4
from sage.functions.log import ln
cn = [(ln(n)**4) for n in range(1, 11)]
cn_conclusion = "The sequence cn is defined as [ln(n)]^4 and does not converge. It is not a monotone sequence."

# Plot the sequences
plt.figure(figsize=(12, 5))
plt.subplot(131)
plt.plot(an, label='an')
plt.title('Sequence an: 8n!')
plt.xlabel('n')
plt.ylabel('Value')
plt.grid()

plt.subplot(132)
plt.plot(bn, label='bn')
plt.title('Sequence bn: n + 1')
plt.xlabel('n')
plt.ylabel('Value')
plt.grid()

plt.subplot(133)
plt.plot(cn, label='cn')
plt.title('Sequence cn: [ln(n)]^4')
plt.xlabel('n')
plt.ylabel('Value')
plt.grid()

plt.tight_layout()
plt.show()

an_conclusion, bn_conclusion, cn_conclusion

# Initialize a list for partial sums
partial_sums = []
sum_terms = 0

# Calculate and store the partial sums
for k in range(1, 21):
    sum_terms += 1/k^2
    partial_sums.append(sum_terms.n())

# Display the partial sums along with their respective term indices
for k, partial_sum in enumerate(partial_sums, start=1):
    print(f"Partial Sum S({k}) = {partial_sum}")

# Create a list of the first 50 terms of the series (-1)^n/n
terms = [(-1)^n/n for n in range(1, 51)]

# Create a list of partial sums up to the 50th term
partial_sums = [sum(terms[:k]).n() for k in range(1, 51)]

# Create a graph to visualize the terms and partial sums
terms_plot = list_plot(terms, color='blue', legend_label='Terms')
partial_sums_plot = list_plot(partial_sums, color='red', legend_label='Partial Sums')

# Display the graph
show(terms_plot + partial_sums_plot)

# Import relevant functions and libraries
from sage.symbolic.integration.integral import definite_integral
from sage.plot.plot import plot

# Define symbolic variables
i, j, k, n = var('i j k n')

# Function to handle divergent cases
def handle_divergent_sum(expression, variable, start, end):
    result = None
    result_numerical = None
    try:
        result = sum(expression, variable, start, end)
        if not result.is_infinite():
            result_numerical = result.n()
    except Exception as e:
        result = "Divergent (Non-convergent Series)"
        result_numerical = "Divergent (Non-convergent Series)"
    return result, result_numerical

# Summation i: Summing numbers with a formula (1/1^2, 1/2^2, 1/3^2, ...)
sum_i, sum_i_numerical = handle_divergent_sum(1/i^2, i, 1, infinity)

# Summation ii: Another series (1/(2^2 - 1), 1/(3^2 - 1), 1/(4^2 - 1), ...)
sum_ii = sum(1/(j^2 - 1), j, 2, infinity)
sum_ii_strategies = "Strategies to calculate this sum may include partial fraction decomposition."

# Summation iii: Similar series with a different formula
sum_iii = sum(1/(k^2 - 3), k, 1, infinity)
sum_iii_strategies = "Strategies to calculate this sum may include partial fraction decomposition."

# Summation iv: Summing numbers with a formula (1/1, 1/2, 1/3, ...)
sum_iv, sum_iv_numerical = handle_divergent_sum(1/n, n, 1, infinity)

# Problematic sum v: A more complex sum with the natural logarithm (ln)
problematic_sum = None
problematic_sum_numerical = None
try:
    problematic_sum = sum(ln((n+1)^2)/n^2, n, 1, infinity)
    if not problematic_sum.is_infinite():
        numerical_result, _ = numerical_integral(problematic_sum, 1, infinity)
        problematic_sum_numerical = numerical_result
except Exception as e:
    problematic_sum_numerical = "positive infinity"

# Calculate an integral
x = var('x')
integral_result = definite_integral(ln((x+1)^2)/x^2, x, 1, infinity)

# Define a function for right Riemann approximation and plot it
def right_riemann_plot(f, a, b, n):
    p = plot(f, a, b, thickness=2)
    deltax = (b - a) / n
    oline = []
    for j in range(n):
        oline += [(a + j * deltax, 0), (a + j * deltax, f(a + (j+1) * deltax)),
                  (a + (j+1) * deltax, f(a + (j+1) * deltax)), (a + (j+1) * deltax, 0)]
    pline = plot(line(oline, rgbcolor='black'))
    show(p + pline)

# Plot a right Riemann approximation for a function
f(x) = (x - 3) * (x - 5) * (x - 7) + 40
right_riemann_plot(f, 0, 3, 4)

# Organize results into a summary with labels and explanations
result_summary = [
    "Result of sum_i: The result of ∑∞ 1/i^2 is", sum_i,
    "Result of sum_ii:", sum_ii_strategies, " The value of ∑∞ 1/((j^2 - 1)) is", sum_ii,
    "Result of sum_iii:", sum_iii_strategies, "The value of ∑∞ 1/(k^2 - 3) is", sum_iii,
    "Result of sum_iv: The result of ∑∞ 1/n is", sum_iv,
    "Result of problematic_sum: Symbolic expression for a problematic sum is", problematic_sum,
    "Numerical approximation of problematic_sum is", problematic_sum_numerical,
    "Result of integral_result: The integral of ln((x+1)^2)/x^2 from 1 to infinity is", integral_result
]

# Filter out None and NaN results
result_summary = [item for item in result_summary if item is not None]

# Display the result summary with labels and explanations
result_summary

# Define the variable and the series
n = var('n')
series = sum(ln((n+1)^2)/n^2, n, 1, 50)

# Calculate the numerical approximation of s50
numerical_approximation = series.n()

# Calculate the remainder using numerical integration
remainder, _ = numerical_integral(ln((n+1)^2)/n^2, 50, infinity)

# Calculate the error interval
error_lower_bound = (numerical_approximation - remainder).n()
error_upper_bound = (numerical_approximation + remainder).n()

# Display the result with an explanation
result = f"The numerical approximation of s50 for the series ∑(ln((n+1)^2)/n^2) from n=1 to 50 is approximately {numerical_approximation}."
result

# Display the error interval
error_interval = f"The error interval is approximately [{error_lower_bound}, {error_upper_bound}]."
error_interval

{
    "Result Explanation": result,
    "Error Interval": error_interval
}

sageStuff='Hi'
sageStuff

# Define the variable and the series
n = var('n')
series = sum((5*n + 1)/(n^3 + n), n, 1, 4)

# Calculate the numerical approximation of s4
s4 = series.n()

# Define the function to graph
f(x) = (5*x + 1)/(x^3 + x)

# Create a symbolic expression for the function
f_symbolic = (5*x + 1)/(x^3 + x)

# Calculate the integral of the function from 1 to infinity
integral_result = definite_integral(f_symbolic, x, 1, infinity)

# Define a function for the right Riemann approximation and plot it
def right_riemann_plot(f, a, b, n):
    p = plot(f, a, b, thickness=2)
    deltax = (b - a) / n
    oline = []
    for j in range(n):
        oline += [(a + j * deltax, 0), (a + j * deltax, f(a + (j+1) * deltax)),
                  (a + (j+1) * deltax, f(a + (j+1) * deltax)), (a + (j+1) * deltax, 0)]
    pline = plot(line(oline, rgbcolor='black'))
    show(p + pline)

# Plot the right Riemann approximation
right_riemann_plot(f, 1, 4, 4)

# Estimate the error
error_lower_bound = (s4 - integral_result).n()
error_upper_bound = (s4 + integral_result).n()

# Display the results with explanations
result_explanation = f"The numerical approximation of s4 for the series is approximately {s4}."
integral_explanation = f"The integral result (from 1 to infinity) is approximately {integral_result}."
error_explanation = f"The error is estimated to be between {error_lower_bound} and {error_upper_bound}. The integral essentially represents the exact accumulation (area under the curve) of the function from 1 to infinity. The difference between the numerical approximation and the integral gives an estimate of the error introduced by simplifying the series at a finite point (in this case, at n=4). By adding and subtracting the integral result, we create a range within which we expect the true sum of the series to lie. This range serves as the upper and lower bounds on the error."

{
    "Result Explanation": result_explanation,
    "Integral Explanation": integral_explanation,
    "Error Explanation": error_explanation
}