diff --git a/homework10/Homework10_Q1.py b/homework10/Homework10_Q1.py new file mode 100644 index 0000000..a73f878 --- /dev/null +++ b/homework10/Homework10_Q1.py @@ -0,0 +1,249 @@ +import numpy as np +import matplotlib.pyplot as plt +import sys + +# Configuration for plotting +plt.style.use('seaborn-v0_8') +plt.rcParams['figure.figsize'] = (18, 6) + +def newton_method(f, df, x0, tolerance=1e-12, max_iter=20): + """ + Implements Newton's Method with overflow protection. + Returns: list of iteration history [(iter, x, error)] + """ + history = [] + x = x0 + + for i in range(max_iter): + try: + # Check for divergence first + if abs(x) > 1e100: + break + + fx = f(x) + dfx = df(x) + + if abs(dfx) < 1e-15: # Avoid division by zero + break + + history.append(x) + + x_new = x - fx / dfx + + # Check convergence + if abs(x_new - x) < tolerance: + x = x_new + history.append(x) + break + + x = x_new + except OverflowError: + # Catch errors if numbers become too large + break + + return history + +def chord_method(f, df_x0, x0, tolerance=1e-12, max_iter=50): + """ + Implements Chord Method using fixed slope m = f'(x0) with overflow protection. + Returns: list of iteration history + """ + history = [] + x = x0 + m = df_x0 # Fixed slope + + for i in range(max_iter): + try: + if abs(x) > 1e100: + break + + fx = f(x) + history.append(x) + + if abs(m) < 1e-15: + break + + x_new = x - fx / m + + if abs(x_new - x) < tolerance: + x = x_new + history.append(x) + break + x = x_new + except OverflowError: + break + + return history + +def secant_method(f, x0, x_minus_1, tolerance=1e-12, max_iter=20): + """ + Implements Secant Method with overflow protection. + Returns: list of iteration history + """ + history = [] + x_curr = x0 + x_prev = x_minus_1 + + for i in range(max_iter): + try: + if abs(x_curr) > 1e100: + break + + history.append(x_curr) + + fx_curr = f(x_curr) + fx_prev = f(x_prev) + + if abs(fx_curr - fx_prev) < 1e-15: + break + + # Secant update + x_new = x_curr - fx_curr * (x_curr - x_prev) / (fx_curr - fx_prev) + + if abs(x_new - x_curr) < tolerance: + x_curr = x_new + history.append(x_curr) + break + + x_prev = x_curr + x_curr = x_new + except OverflowError: + break + + return history + +def solve_and_save(): + # Definition of problems + problems = [ + { + 'id': 'a', + 'title': r'$f(x) = \cos(x) - x, x_0 = 0.5$', + 'func_name': 'cos(x) - x', + 'f': lambda x: np.cos(x) - x, + 'df': lambda x: -np.sin(x) - 1, + 'x0': 0.5, + 'true_root': 0.73908513321516064165531208767387340401341175890075746496568063577328465488354759 # From wolframalpha + }, + { + 'id': 'b', + 'title': r'$f(x) = \sin(x), x_0 = 3$', + 'func_name': 'sin(x)', + 'f': lambda x: np.sin(x), + 'df': lambda x: np.cos(x), + 'x0': 3.0, + 'true_root': np.pi + }, + { + 'id': 'c', + 'title': r'$f(x) = x^2 + 1, x_0 = 0.5$', + 'func_name': 'x^2 + 1', + 'f': lambda x: x**2 + 1, + 'df': lambda x: 2*x, + 'x0': 0.5, + 'true_root': None # No real root + } + ] + + # Open file for writing text results + output_filename = 'numerical_results.txt' + image_filename = 'convergence_results.png' + + with open(output_filename, 'w', encoding='utf-8') as txt_file: + + fig, axes = plt.subplots(1, 3) + + for idx, prob in enumerate(problems): + header_str = f"\n{'='*60}\nProblem ({prob['id']}): {prob['func_name']}, x0={prob['x0']}\n{'='*60}\n" + print(header_str) # Print to console to show progress + txt_file.write(header_str) + + x0 = prob['x0'] + f = prob['f'] + df = prob['df'] + + # 1. Run Newton + hist_newton = newton_method(f, df, x0) + + # 2. Run Chord (m = f'(x0)) + hist_chord = chord_method(f, df(x0), x0) + + # 3. Run Secant (x_-1 = 0.99 * x0) + x_minus_1 = 0.99 * x0 + hist_secant = secant_method(f, x0, x_minus_1) + + methods_data = { + "Newton": hist_newton, + "Chord": hist_chord, + "Secant": hist_secant + } + + ax = axes[idx] + + for m_name, history in methods_data.items(): + iters = list(range(len(history))) + errors = [] + + # Calculate errors for plotting + # Note: Handling overflow for plotting large errors in problem c + for x in history: + try: + if prob['true_root'] is not None: + val = abs(x - prob['true_root']) + else: + val = abs(f(x)) # Residual for problem c + errors.append(val) + except OverflowError: + errors.append(float('inf')) + + # Write to TXT file + txt_file.write(f"\nMethod: {m_name}\n") + txt_file.write(f"{'Iter':<5} | {'Value (x_k)':<25} | {'Error/Resid':<25}\n") + txt_file.write("-" * 60 + "\n") + + for i, val in enumerate(history): + # Only write full history if short, or head/tail if long + if len(history) < 20 or i < 5 or i > len(history) - 3: + if prob['true_root'] is not None: + err_val = abs(val - prob['true_root']) + else: + try: + err_val = abs(f(val)) + except OverflowError: + err_val = float('inf') + + txt_file.write(f"{i:<5} | {val:<25.10g} | {err_val:<25.6e}\n") + elif i == 5: + txt_file.write("...\n") + + # Add to plot + if prob['true_root'] is not None: + ylabel_text = "Error |x_k - x*|" + log_scale = True + ax.plot(iters, errors, marker='o', label=m_name) + else: + ylabel_text = "Residual |f(x_k)|" + log_scale = False + # For problem C, data explodes, so we limit plotting range to avoid ruining the graph + valid_errors = [e for e in errors if e < 1e10] + valid_iters = iters[:len(valid_errors)] + ax.plot(valid_iters, valid_errors, marker='x', label=m_name) + + ax.set_title(f"Problem ({prob['id']})\n{prob['title']}") + ax.set_xlabel("Iteration") + ax.set_ylabel(ylabel_text) + if log_scale: + ax.set_yscale('log') + ax.legend() + if prob['id'] == 'c': + ax.set_ylim(0, 20) + ax.set_xticks(range(0, len(iters), 1)) + ax.grid(True) + + print(f"\nText results saved to: {output_filename}") + + plt.tight_layout() + plt.savefig(image_filename) + print(f"Chart saved to: {image_filename}") + +if __name__ == "__main__": + solve_and_save() diff --git a/homework10/comment.txt b/homework10/comment.txt new file mode 100644 index 0000000..59112df --- /dev/null +++ b/homework10/comment.txt @@ -0,0 +1 @@ +For problems (a) and (b), which have real roots, all methods converged; Newton's method was the fastest (quadratic), followed by the Secant method (superlinear), with the Chord method being the slowest (linear). However, for problem (c) (x^2 + 1), which has no real roots, all methods failed: the Chord method rapidly diverged to infinity, while Secant and Newton's method oscillated. \ No newline at end of file diff --git a/homework10/convergence_results.png b/homework10/convergence_results.png new file mode 100644 index 0000000..e906462 Binary files /dev/null and b/homework10/convergence_results.png differ diff --git a/homework10/numerical_results.txt b/homework10/numerical_results.txt new file mode 100644 index 0000000..43691c0 --- /dev/null +++ b/homework10/numerical_results.txt @@ -0,0 +1,112 @@ + +============================================================ +Problem (a): cos(x) - x, x0=0.5 +============================================================ + +Method: Newton +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 2.390851e-01 +1 | 0.7552224171 | 1.613728e-02 +2 | 0.7391416661 | 5.653293e-05 +3 | 0.7390851339 | 7.056461e-10 +4 | 0.7390851332 | 0.000000e+00 +5 | 0.7390851332 | 0.000000e+00 + +Method: Chord +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 2.390851e-01 +1 | 0.7552224171 | 1.613728e-02 +2 | 0.7369022576 | 2.182876e-03 +3 | 0.7393704622 | 2.853290e-04 +4 | 0.7390476612 | 3.747205e-05 +... +13 | 0.7390851332 | 4.333200e-13 +14 | 0.7390851332 | 5.695444e-14 + +Method: Secant +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 2.390851e-01 +1 | 0.7556018135 | 1.651668e-02 +2 | 0.738106944 | 9.781892e-04 +3 | 0.7390815948 | 3.538433e-06 +4 | 0.739085134 | 7.646579e-10 +5 | 0.7390851332 | 5.551115e-16 +6 | 0.7390851332 | 0.000000e+00 + +============================================================ +Problem (b): sin(x), x0=3.0 +============================================================ + +Method: Newton +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 3 | 1.415927e-01 +1 | 3.142546543 | 9.538895e-04 +2 | 3.141592653 | 2.893161e-10 +3 | 3.141592654 | 0.000000e+00 +4 | 3.141592654 | 0.000000e+00 + +Method: Chord +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 3 | 1.415927e-01 +1 | 3.142546543 | 9.538895e-04 +2 | 3.141583011 | 9.642404e-06 +3 | 3.141592751 | 9.747184e-08 +4 | 3.141592653 | 9.853101e-10 +5 | 3.141592654 | 9.960477e-12 +6 | 3.141592654 | 1.003642e-13 +7 | 3.141592654 | 1.332268e-15 + +Method: Secant +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 3 | 1.415927e-01 +1 | 3.142873442 | 1.280788e-03 +2 | 3.141588403 | 4.250744e-06 +3 | 3.141592654 | 1.158629e-12 +4 | 3.141592654 | 0.000000e+00 +5 | 3.141592654 | 0.000000e+00 + +============================================================ +Problem (c): x^2 + 1, x0=0.5 +============================================================ + +Method: Newton +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 1.250000e+00 +1 | -0.75 | 1.562500e+00 +2 | 0.2916666667 | 1.085069e+00 +3 | -1.568452381 | 3.460043e+00 +4 | -0.4654406117 | 1.216635e+00 +... +18 | -1.820817243 | 4.315375e+00 +19 | -0.6358066524 | 1.404250e+00 + +Method: Chord +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 1.250000e+00 +1 | -0.75 | 1.562500e+00 +2 | -2.3125 | 6.347656e+00 +3 | -8.66015625 | 7.599831e+01 +4 | -84.65846252 | 7.168055e+03 +... +8 | -7.660265913e+30 | 5.867967e+61 +9 | -5.867967385e+61 | 3.443304e+123 + +Method: Secant +Iter | Value (x_k) | Error/Resid +------------------------------------------------------------ +0 | 0.5 | 1.250000e+00 +1 | -0.756281407 | 1.571962e+00 +2 | 5.37745098 | 2.991698e+01 +3 | -1.096446714 | 2.202195e+00 +4 | -1.610857647 | 3.594862e+00 +... +18 | 0.3478146096 | 1.120975e+00 +19 | -2.615940634 | 7.843145e+00