araff/NumericalMethodsSolver

1 2import gradio as gr3import numpy as np4import pandas as pd5import matplotlib.pyplot as plt6import sympy as sp7 8# ============================================================9# Utility functions10# ============================================================11 12def make_rhs_function(rhs_str, order):13    x = sp.symbols('x')14    y_symbols = sp.symbols(f'y0:{order}')15    expr = sp.sympify(rhs_str)16    return sp.lambdify((x, *y_symbols), expr, modules=['numpy'])17 18def make_exact_function(expr_str):19    x = sp.symbols('x')20    expr = sp.sympify(expr_str)21    return sp.lambdify(x, expr, modules=['numpy'])22 23def system_rhs(x, Y, g, order):24    F = np.zeros(order, dtype=float)25    for i in range(order - 1):26        F[i] = Y[i + 1]27    F[order - 1] = g(x, *Y)28    return F29 30# ============================================================31# Numerical methods32# ============================================================33 34def forward_euler(g, order, x0, Y0, x_end, h):35    nsteps = int(round((x_end - x0) / h))36    x = np.linspace(x0, x_end, nsteps + 1)37    Y = np.zeros((nsteps + 1, order), dtype=float)38    Y[0, :] = Y039 40    for i in range(nsteps):41        Y[i + 1, :] = Y[i, :] + h * system_rhs(x[i], Y[i, :], g, order)42 43    return x, Y44 45def modified_euler(g, order, x0, Y0, x_end, h):46    nsteps = int(round((x_end - x0) / h))47    x = np.linspace(x0, x_end, nsteps + 1)48    Y = np.zeros((nsteps + 1, order), dtype=float)49    Y[0, :] = Y050 51    for i in range(nsteps):52        k1 = system_rhs(x[i], Y[i, :], g, order)53        Y_pred = Y[i, :] + h * k154        k2 = system_rhs(x[i + 1], Y_pred, g, order)55        Y[i + 1, :] = Y[i, :] + (h / 2.0) * (k1 + k2)56 57    return x, Y58 59def rk4(g, order, x0, Y0, x_end, h):60    nsteps = int(round((x_end - x0) / h))61    x = np.linspace(x0, x_end, nsteps + 1)62    Y = np.zeros((nsteps + 1, order), dtype=float)63    Y[0, :] = Y064 65    for i in range(nsteps):66        k1 = system_rhs(x[i], Y[i, :], g, order)67        k2 = system_rhs(x[i] + h / 2.0, Y[i, :] + h * k1 / 2.0, g, order)68        k3 = system_rhs(x[i] + h / 2.0, Y[i, :] + h * k2 / 2.0, g, order)69        k4 = system_rhs(x[i] + h, Y[i, :] + h * k3, g, order)70 71        Y[i + 1, :] = Y[i, :] + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)72 73    return x, Y74 75def backward_euler(g, order, x0, Y0, x_end, h, tol=1e-10, max_iter=100):76    nsteps = int(round((x_end - x0) / h))77    x = np.linspace(x0, x_end, nsteps + 1)78    Y = np.zeros((nsteps + 1, order), dtype=float)79    Y[0, :] = Y080 81    for i in range(nsteps):82        x_next = x[i + 1]83        guess = Y[i, :] + h * system_rhs(x[i], Y[i, :], g, order)84 85        for _ in range(max_iter):86            new_guess = Y[i, :] + h * system_rhs(x_next, guess, g, order)87            if np.linalg.norm(new_guess - guess, ord=np.inf) < tol:88                guess = new_guess89                break90            guess = new_guess91 92        Y[i + 1, :] = guess93 94    return x, Y95 96# ============================================================97# Main solver98# ============================================================99 100def solve_ode_problem(order, rhs_expr, ic_text, x0, x_end, h, methods, exact_expr=''):101    order = int(order)102 103    if order < 1:104        raise ValueError('Order must be at least 1.')105 106    if h <= 0:107        raise ValueError('Step size h must be positive.')108 109    if x_end <= x0:110        raise ValueError('x_end must be greater than x0.')111 112    ratio = (x_end - x0) / h113    if abs(ratio - round(ratio)) > 1e-10:114        raise ValueError('Choose h so that (x_end - x0)/h is an integer.')115 116    ic_values = [float(v.strip()) for v in ic_text.split(',') if v.strip() != '']117    if len(ic_values) != order:118        raise ValueError(119            f'You entered {len(ic_values)} initial condition(s), but order = {order}. '120            f'Please provide exactly {order} initial conditions.'121        )122 123    if len(methods) == 0:124        raise ValueError('Select at least one numerical method.')125 126    g = make_rhs_function(rhs_expr, order)127    Y0 = np.array(ic_values, dtype=float)128 129    results = {}130 131    if 'Forward Euler' in methods:132        results['Forward Euler'] = forward_euler(g, order, x0, Y0, x_end, h)133 134    if 'Backward Euler' in methods:135        results['Backward Euler'] = backward_euler(g, order, x0, Y0, x_end, h)136 137    if 'Modified Euler (Heun)' in methods:138        results['Modified Euler (Heun)'] = modified_euler(g, order, x0, Y0, x_end, h)139 140    if 'RK4' in methods:141        results['RK4'] = rk4(g, order, x0, Y0, x_end, h)142 143    x_vals = list(results.values())[0][0]144    df = pd.DataFrame({'x': x_vals})145 146    for method_name, (_, Y) in results.items():147        for j in range(order):148            if j == 0:149                col_name = f'{method_name} : y'150            elif j == 1:151                col_name = f"{method_name} : y'"152            elif j == 2:153                col_name = f"{method_name} : y''"154            else:155                col_name = f'{method_name} : d^{j}y/dx^{j}'156            df[col_name] = Y[:, j]157 158    exact_used = False159    if exact_expr.strip():160        exact_f = make_exact_function(exact_expr)161        y_exact = exact_f(x_vals)162        df['Exact y'] = y_exact163        exact_used = True164 165        for method_name, (_, Y) in results.items():166            df[f'Error ({method_name})'] = np.abs(Y[:, 0] - y_exact)167 168    csv_path = 'solution_results.csv'169    df.to_csv(csv_path, index=False)170 171    fig1, ax1 = plt.subplots(figsize=(9, 5))172    for method_name, (_, Y) in results.items():173        ax1.plot(x_vals, Y[:, 0], marker='o', label=method_name)174 175    if exact_used:176        ax1.plot(x_vals, df['Exact y'], linestyle='--', linewidth=2, label='Exact y')177 178    ax1.set_xlabel('x')179    ax1.set_ylabel('y')180    ax1.set_title('Comparison of Numerical Methods')181    ax1.grid(True)182    ax1.legend()183 184    fig2 = None185    if order > 1:186        fig2, ax2 = plt.subplots(figsize=(9, 5))187        for method_name, (_, Y) in results.items():188            for j in range(1, order):189                if j == 1:190                    lab = f"{method_name} : y'"191                elif j == 2:192                    lab = f"{method_name} : y''"193                else:194                    lab = f'{method_name} : order {j}'195                ax2.plot(x_vals, Y[:, j], marker='o', label=lab)196 197        ax2.set_xlabel('x')198        ax2.set_ylabel('Derivative values')199        ax2.set_title('Higher-Order State Variables')200        ax2.grid(True)201        ax2.legend()202 203    summary = (204        f'Order: {order}\n'205        f'Equation entered: y^({order}) = {rhs_expr}\n'206        f'Initial conditions at x0 = {x0}: {ic_text}\n'207        f'x_end = {x_end}\n'208        f'Step size h = {h}\n'209        f'Methods used: {', '.join(methods)}\n'210        f'Exact solution provided: {'Yes' if exact_used else 'No'}'211    )212 213    return df, fig1, fig2, summary, csv_path214 215DESCRIPTION = '''216### Description217This app solves ordinary differential equations written in the form:218 219**y^(n) = f(x, y, y\', y\'', ..., y^(n-1))**220 221### How to enter the problem222- Enter the **order** of the differential equation.223- Enter the **right-hand side** only.224- Use:225  - `y0` for `y`226  - `y1` for `y'`227  - `y2` for `y''`228  - and so on229- Enter initial conditions as comma-separated values at `x = x0`.230 231For an equation of order `n`, enter exactly `n` initial conditions:232- `y(x0), y'(x0), y''(x0), ..., y^(n-1)(x0)`233 234### Examples235**First-order**236- Equation: y' = y + 2x - x^2237- Order: `1`238- RHS: `y0 + 2*x - x**2`239- Initial conditions: `1`240 241**Second-order**242- Equation: y'' = -y243- Order: `2`244- RHS: `-y0`245- Initial conditions: `0,1`246 247**Third-order**248- Equation: y^(3) = x + y2 - 2*y1 + y0249- Order: `3`250- RHS: `x + y2 - 2*y1 + y0`251- Initial conditions: `1,0,-1`252 253### Methods available254- Forward Euler255- Backward Euler256- Modified Euler (Heun)257- Runge-Kutta 4th Order (RK4)258 259Developed by **Jude Okolie**260'''261 262with gr.Blocks() as demo:263    gr.Markdown('# Numerical Methods Solver for Differential Equations')264    gr.Markdown(DESCRIPTION)265 266    with gr.Row():267        order = gr.Number(label='Order of ODE', value=1, precision=0)268        x0 = gr.Number(label='x0', value=0.0)269        x_end = gr.Number(label='x_end', value=1.5)270        h = gr.Number(label='Step size h', value=0.375)271 272    rhs_expr = gr.Textbox(273        label='Right-hand side f(x, y, y\', ...)',274        value='y0 + 2*x - x**2',275        lines=1276    )277 278    ic_text = gr.Textbox(279        label='Initial conditions at x = x0 (comma-separated)',280        value='1',281        lines=1282    )283 284    exact_expr = gr.Textbox(285        label='Optional exact solution y(x)',286        value='x**2 + exp(x)',287        lines=1288    )289 290    methods = gr.CheckboxGroup(291        choices=['Forward Euler', 'Backward Euler', 'Modified Euler (Heun)', 'RK4'],292        value=['Forward Euler', 'Backward Euler', 'Modified Euler (Heun)', 'RK4'],293        label='Select numerical methods'294    )295 296    solve_btn = gr.Button('Solve')297 298    results_table = gr.Dataframe(label='Numerical Results')299    plot_y = gr.Plot(label='Solution Plot')300    plot_derivs = gr.Plot(label='Higher-Order State Variables')301    summary = gr.Textbox(label='Solution Summary', lines=8)302    csv_file = gr.File(label='Download Results CSV')303 304    solve_btn.click(305        fn=solve_ode_problem,306        inputs=[order, rhs_expr, ic_text, x0, x_end, h, methods, exact_expr],307        outputs=[results_table, plot_y, plot_derivs, summary, csv_file]308    )309 310if __name__ == '__main__':311    demo.launch()312