← All virtual labs
Numerical Methods · Virtual lab

Roots, interpolation, quadrature, ODEs & linear systems

Six self-contained experiments, each a genuinely working numerical engine written in plain JavaScript — no libraries. You type a real function f(x) (parsed by a hand-rolled tokenizer and recursive-descent evaluator) and every control drives the computation: iterations animate, convergence and error are tabulated per step, and the exact value is compared against the numerical estimate. Pick an experiment on the left, then move through Aim, Theory, Procedure, the live Simulation, a graded Self-assessment, and References.

1 · Root finding

Bisection, Newton-Raphson, secant and false-position on your own f(x), with animated iterates and per-step error
To solve f(x) = 0 for a chosen continuous function using four classical iterative methods, to watch each iterate converge on the canvas, and to compare their convergence orders by tabulating the approximate error at every step.

Bisection brackets a root between a and b with f(a)·f(b) < 0 and halves the interval each step, keeping the half that still changes sign. It always converges (linearly), gaining one bit per iteration: the error bound is (b−a)/2 to the power n.

Bisection: c = (a + b) / 2, keep [a,c] if f(a)f(c) < 0 else [c,b]

Newton-Raphson uses the tangent line. From the Taylor expansion, the next iterate is the x-intercept of the tangent at the current point. It converges quadratically near a simple root but needs the derivative and can diverge.

Newton: x(n+1) = x(n) − f(x(n)) / f'(x(n))

The secant method replaces f' by a finite-difference slope through the last two points, giving superlinear order about 1.618 without needing a derivative. False position (regula falsi) uses the same secant formula but, like bisection, keeps a sign-changing bracket, so it is guaranteed to converge.

Secant / false-position: x(n+1) = x(n) − f(x(n))·(x(n) − x(n−1)) / (f(x(n)) − f(x(n−1)))

The approximate relative error reported each step is |x(n+1) − x(n)| / |x(n+1)|; iteration stops when it falls below the tolerance or the iteration cap is reached.

  1. Open Simulation. A default function and bracket are loaded.
  2. Edit f(x) (for example x^3 - x - 2, cos(x) - x, exp(-x) - x). The plot redraws live; a red box means the expression failed to parse.
  3. Choose a method. For bracketing methods set a and b so the curve crosses between them; for open methods set the start point(s).
  4. Set the tolerance and max iterations, then press Run. Watch the iterate markers march toward the root.
  5. Read the per-step table: the current estimate, f(x) and the approximate error. Note how Newton roughly doubles the number of correct digits each step while bisection halves the bracket.
  6. Press Step to advance one iteration at a time, or Reset to clear.

Function & iteratesready

Tip: use ^ for powers, and functions sin cos tan exp log sqrt abs. Constants pi and e are available. The dashed grey line is y = 0.

Controls

Root est.
--
f(root)
--
Iterations
0
Rel. error
--
nabxf(x)err
  • Burden & Faires — Numerical Analysis, 9th ed., Ch. 2 (Solutions of Equations in One Variable). Brooks/Cole.
  • Atkinson — An Introduction to Numerical Analysis, 2nd ed., Sec. 2. Wiley.
  • Virtual Labs (IIT) — Numerical Analysis: Root Finding, nalwss-coep.vlabs.ac.in.

2 · Polynomial interpolation

Lagrange and Newton divided-difference polynomials through your data points, drawn as a continuous curve
To construct the unique polynomial of degree at most n through n+1 given points using both the Lagrange and the Newton divided-difference forms, to confirm the two forms give the same curve, and to estimate the value of the interpolant at any query point.

Given n+1 distinct points there is exactly one polynomial of degree at most n passing through all of them. The Lagrange form writes it as a weighted sum of basis polynomials, each of which is one at its own node and zero at every other node.

P(x) = sum over i of y(i) · L(i, x), L(i, x) = product over j not equal to i of (x − x(j)) / (x(i) − x(j))

The Newton divided-difference form builds the same polynomial incrementally, so adding a point costs only one extra term. The coefficients are the leading divided differences computed from a triangular table.

P(x) = f[x0] + f[x0,x1](x − x0) + f[x0,x1,x2](x − x0)(x − x1) + ...
f[x(i),...,x(i+k)] = ( f[x(i+1),...] − f[x(i),...] ) / ( x(i+k) − x(i) )

Both are exact at the nodes; between them they can oscillate (Runge phenomenon) for high degree on equally spaced data. The interpolation error involves the (n+1)-th derivative of the true function and the node product.

  1. Open Simulation. Several points are preloaded and the interpolating polynomial is drawn through them.
  2. Click on the canvas to add a point, or use x, y and Add point. Use Random or Clear to reshape the data set.
  3. Toggle Lagrange and Newton; both curves should coincide exactly — they are the same polynomial in two algebraic dresses.
  4. Enter a query x and read the interpolated value; a vertical marker shows it on the curve.
  5. Read the printed divided-difference table — its top row gives the Newton coefficients.

Interpolant0 points

Click the plot to drop a point. The polynomial has degree = points − 1. Watch for wild swings between widely spaced nodes — that is the Runge phenomenon.

Data & query

Both
Lagrange
Newton
P = --
Points
0
Degree
--
Divided-difference table (Newton coefficients along the top diagonal):
  • Burden & Faires — Numerical Analysis, 9th ed., Ch. 3 (Interpolation and Polynomial Approximation). Brooks/Cole.
  • Hildebrand — Introduction to Numerical Analysis, 2nd ed., Ch. 2. Dover.
  • Virtual Labs (IIT) — Numerical Analysis: Interpolation, nalwss-coep.vlabs.ac.in.

3 · Numerical integration

Trapezoidal, Simpson 1/3 and Simpson 3/8 quadrature with drawn panels and error against the reference value
To approximate a definite integral of f(x) over [a,b] by the composite trapezoidal, Simpson one-third and Simpson three-eighths rules, to visualise the approximating panels, and to study how the error shrinks as the number of subintervals grows.

A Newton-Cotes rule replaces f by an interpolating polynomial on each panel and integrates that exactly. With h = (b−a)/n the composite trapezoidal rule joins consecutive points by straight lines:

T = (h/2) [ f0 + 2(f1 + f2 + ... + f(n−1)) + fn ], error ~ O(h squared)

Simpson's 1/3 rule fits a parabola over each pair of panels (so n must be even) and is exact for cubics:

S1/3 = (h/3) [ f0 + 4(f1+f3+...) + 2(f2+f4+...) + fn ], error ~ O(h to the 4)

Simpson's 3/8 rule fits a cubic over each group of three panels (n a multiple of 3):

S3/8 = (3h/8) [ f0 + 3(f1+f2) + 3(f4+f5) + ... + 2(f3+f6+...) + fn ]

The reported error is the absolute difference from a high-resolution reference value (Simpson with many points), so you can watch the trapezoidal error fall like h squared while Simpson falls like h to the fourth.

  1. Open Simulation. The default integrand and limits are loaded and the panels are shaded.
  2. Edit f(x), the limits a and b, and the number of subintervals n with the slider.
  3. Pick a rule. For Simpson 1/3, n is rounded up to an even number; for 3/8, to a multiple of three (the panel says which).
  4. Read the estimate, the reference value and the absolute error.
  5. Increase n and watch the error drop — the panel reports the empirical order of accuracy as you do.

Quadrature panelsready

Shaded regions are the panels actually summed. Straight tops mean trapezoidal; curved tops mean the fitted parabola or cubic of Simpson's rules.

Controls

8
Trapz
Simp 1/3
Simp 3/8
Estimate
--
Reference
--
Abs. error
--
Step h
--
Effective n: -- · empirical order p: --
  • Burden & Faires — Numerical Analysis, 9th ed., Ch. 4 (Numerical Differentiation and Integration). Brooks/Cole.
  • Davis & Rabinowitz — Methods of Numerical Integration, 2nd ed. Academic Press.
  • Virtual Labs (IIT) — Numerical Analysis: Numerical Integration, nalwss-coep.vlabs.ac.in.

4 · ODE solvers

Euler, Heun (improved Euler) and classical RK4 on dy/dx = f(x,y), against the analytic solution
To integrate a first-order initial-value problem by the explicit Euler, Heun and fourth-order Runge-Kutta methods, to plot each numerical solution against the known analytic curve, and to compare their global errors as the step size changes.

For dy/dx = f(x,y) with y(x0) = y0, a one-step method advances by h. Euler's method follows the slope at the start of the step — it is first order, so halving h roughly halves the error.

Euler: y(n+1) = y(n) + h · f(x(n), y(n))

Heun's method (improved Euler / RK2) averages the slope at the start with a predicted slope at the end, giving second order:

k1 = f(x, y), k2 = f(x+h, y + h·k1)
y(n+1) = y(n) + (h/2)(k1 + k2)

Classical RK4 takes a weighted average of four slope samples and is fourth order — halving h cuts the error by sixteen:

k1 = f(x, y), k2 = f(x+h/2, y+h·k1/2)
k3 = f(x+h/2, y+h·k2/2), k4 = f(x+h, y+h·k3)
y(n+1) = y(n) + (h/6)(k1 + 2k2 + 2k3 + k4)

The panel reports the maximum absolute error against the analytic solution; choose a preset problem so the exact curve is known.

  1. Open Simulation and pick a preset ODE — each comes with its analytic solution.
  2. Set the initial point x0, y0, the end x and the step h (or the number of steps).
  3. Enable the methods you want to compare. The analytic curve is the dashed green line; numerical solutions are coloured polylines.
  4. Read each method's max error. Halve h and confirm Euler error halves, Heun quarters, RK4 drops about sixteen-fold.
  5. Try a larger h to see Euler drift away while RK4 still tracks the true curve.

Solution curvesready

— — analytic   — Euler   — Heun   — RK4

Controls

Euler
Heun
RK4
Steps
--
y(end) RK4
--
Euler max err
--
Heun max err
--
RK4 max err
--
Exact y(end)
--
  • Burden & Faires — Numerical Analysis, 9th ed., Ch. 5 (Initial-Value Problems for ODEs). Brooks/Cole.
  • Hairer, Norsett & Wanner — Solving Ordinary Differential Equations I, 2nd ed. Springer.
  • Virtual Labs (IIT) — Numerical Analysis: ODE Solvers, nalwss-coep.vlabs.ac.in.

5 · Linear systems

Gauss elimination with partial pivoting, LU decomposition and Gauss-Seidel iteration, with the steps shown
To solve a square linear system Ax = b by Gaussian elimination with partial pivoting and by LU factorisation, then to solve it iteratively with Gauss-Seidel, displaying the intermediate matrices and the iterates so the mechanics of each method are visible.

Gaussian elimination uses row operations to reduce A to upper-triangular form, then back-substitutes. Partial pivoting swaps in the row with the largest pivot magnitude at each column to limit round-off growth and avoid division by a tiny pivot.

multiplier m = a(i,k) / a(k,k); row i := row i − m · row k

LU decomposition records those multipliers in a unit lower-triangular L and the result in an upper-triangular U so that PA = LU. Solving then needs only a forward solve Ly = Pb and a back solve Ux = y — cheap to repeat for many right-hand sides.

A = L U; solve L y = b (forward), then U x = y (backward)

Gauss-Seidel is iterative: it updates each unknown in place using the most recent values of the others. It converges when A is diagonally dominant; the spectral radius of the iteration matrix governs the rate.

x(i) := ( b(i) − sum over j≠i of a(i,j) x(j) ) / a(i,i)

The residual norm ||Ax − b|| measures how well the current solution satisfies the system.

  1. Open Simulation. A 3×3 system is preloaded; edit any entry of A or b.
  2. Choose size (2, 3 or 4) and use Random for a fresh diagonally dominant system, or Reset.
  3. Press Gauss to run elimination with partial pivoting — each pivot and row operation is logged and the triangular matrix is shown.
  4. Press LU to print L and U and the permutation; verify the solution matches Gauss.
  5. Press Gauss-Seidel to iterate; watch the iterates and the residual fall (the badge warns if the matrix is not diagonally dominant).

Working & stepsready

2×2
3×3
4×4

Solution

Method
--
det(A)
--
solution vector x will appear here
Residual ||Ax-b||
--
Diag. dominant
--
  • Golub & Van Loan — Matrix Computations, 4th ed., Ch. 3 (General Linear Systems). Johns Hopkins.
  • Burden & Faires — Numerical Analysis, 9th ed., Ch. 6 & 7. Brooks/Cole.
  • Virtual Labs (IIT) — Numerical Analysis: Systems of Linear Equations, nalwss-coep.vlabs.ac.in.

6 · Least-squares curve fitting

Fit a line or a polynomial of chosen degree to scattered points by solving the normal equations
To fit a polynomial of chosen degree to scattered, noisy data in the least-squares sense by forming and solving the normal equations, and to read off the residual sum of squares and the coefficient of determination R squared.

Least squares chooses coefficients that minimise the sum of squared vertical residuals between the model and the data. For a degree-m polynomial the model is a linear combination of the powers of x, so minimising leads to the normal equations — a square system in the unknown coefficients.

minimise S = sum of ( y(i) − P(x(i)) ) squared
normal equations: (X transpose X) c = X transpose y

Here X is the Vandermonde-style design matrix whose columns are 1, x, x squared, and so on. For a straight line this reduces to the familiar slope and intercept formulas.

line: slope = ( n·Sxy − Sx·Sy ) / ( n·Sxx − Sx squared ), intercept = ( Sy − slope·Sx ) / n

Goodness of fit is summarised by R squared = 1 − SS_res / SS_tot, which is one for a perfect fit and falls toward zero as the model explains less of the variance. Raising the degree always lowers SS_res but risks over-fitting the noise.

  1. Open Simulation. A scatter of noisy points is preloaded and a fit is drawn through it.
  2. Click to add points, or press Generate to sample a known curve plus noise; use Clear to start over.
  3. Set the degree (1 = straight line). The fitted curve and its equation update immediately.
  4. Read the printed coefficients, the residual sum of squares and R squared.
  5. Raise the degree and watch R squared climb toward one while the curve starts to chase individual points — the signature of over-fitting.

Scatter & fit0 points

Click to drop a data point. The purple curve is the least-squares fit; the faint vertical sticks are the residuals being minimised.

Fit controls

1
0.25
Points
0
R squared
--
SS residual
--
RMS error
--
Fitted equation
add at least 2 points
  • Lawson & Hanson — Solving Least Squares Problems. SIAM Classics.
  • Burden & Faires — Numerical Analysis, 9th ed., Sec. 8.1 (Discrete Least Squares Approximation). Brooks/Cole.
  • Virtual Labs (IIT) — Numerical Analysis: Curve Fitting, nalwss-coep.vlabs.ac.in.