Turn math into working code. Implement the algorithms you'll encounter most often in scientific and engineering computation.
#include <stdio.h> #include <math.h> int main(void) { double data[] = {72.5, 85.3, 91.0, 68.2, 77.8, 83.1, 90.5, 65.9}; int n = 8; double sum = 0, sum2 = 0; for (int i = 0; i < n; i++) { sum += data[i]; // sum sum2 += data[i] * data[i]; // sum of squares } double avg = sum / n; double var = sum2 / n - avg * avg; // variance = E[X^2] - (E[X])^2 double sd = sqrt(var); printf("N: %d\n", n); printf("Mean: %.2f\n", avg); printf("Var: %.2f\n", var); printf("Std: %.2f\n", sd); return 0; }
#include <stdio.h> #include <math.h> double f(double x) { return x*x - 2; } // f(x) = x^2 - 2 double df(double x) { return 2*x; } // f'(x) = 2x int main(void) { double x = 2.0; // initial guess double EPS = 1e-10; for (int i = 0; i < 100; i++) { double x_new = x - f(x) / df(x); printf("iter %2d: x = %.15f err = %e\n", i+1, x_new, fabs(x_new - x)); if (fabs(x_new - x) < EPS) break; x = x_new; } printf("sqrt(2) = %.15f\n", x); return 0; }
#include <stdio.h> #include <stdlib.h> #include <time.h> int main(void) { int N; printf("Trials: "); scanf("%d", &N); srand((unsigned)time(NULL)); int inside = 0; for (int i = 0; i < N; i++) { double x = (double)rand() / RAND_MAX; // random in [0,1] double y = (double)rand() / RAND_MAX; if (x*x + y*y <= 1.0) inside++; // inside circle? } double pi = 4.0 * inside / N; printf("pi ≈ %.6f (N=%d)\n", pi, N); printf("true: 3.141593\n"); return 0; }
| Trials N | Estimate (ex.) | Error |
|---|---|---|
| 1,000 | 3.156 | ~0.01 |
| 10,000 | 3.1432 | ~0.001 |
| 100,000 | 3.14052 | ~0.001 |
| 1,000,000 | 3.141284 | ~0.0003 |
| 10,000,000 | 3.141552 | ~0.00004 |
#include <stdio.h> #include <math.h> double f(double x) { return sin(x); } // integrand // Trapezoid rule double trapezoid(double a, double b, int n) { double h = (b - a) / n; double sum = (f(a) + f(b)) / 2.0; for (int i = 1; i < n; i++) sum += f(a + i * h); return h * sum; } // Simpson's rule (higher accuracy) double simpson(double a, double b, int n) { double h = (b - a) / n; double sum = f(a) + f(b); for (int i = 1; i < n; i++) sum += (i % 2 == 0) ? 2*f(a+i*h) : 4*f(a+i*h); return h / 3.0 * sum; } int main(void) { printf("True value of integral of sin(x) from 0 to pi = 2.0\n\n"); int ns[] = {10, 100, 1000}; for (int i = 0; i < 3; i++) { printf("n=%4d trap=%.10f Simpson=%.10f\n", ns[i], trapezoid(0, M_PI, ns[i]), simpson(0, M_PI, ns[i])); } return 0; }
#include <stdio.h> #include <math.h> #define N 3 void gauss(double a[N][N+1], double x[N]) { // forward elimination for (int i = 0; i < N; i++) { for (int j = i+1; j < N; j++) { double ratio = a[j][i] / a[i][i]; for (int k = i; k <= N; k++) a[j][k] -= ratio * a[i][k]; } } // back substitution for (int i = N-1; i >= 0; i--) { x[i] = a[i][N]; for (int j = i+1; j < N; j++) x[i] -= a[i][j] * x[j]; x[i] /= a[i][i]; } } int main(void) { double a[N][N+1] = { { 2, 1, -1, 8}, {-3, -1, 2, -11}, {-2, 1, 2, -3} }; double x[N]; gauss(a, x); for (int i = 0; i < N; i++) printf("x[%d] = %.4f\n", i, x[i]); return 0; // answer: x=2, y=3, z=-1 }
#include <stdio.h> #include <math.h> double f(double x) { return x*x*x - x - 1; } // x^3 - x - 1 = 0 double bisection(double a, double b, double eps) { while (b - a > eps) { double mid = (a + b) / 2.0; if (f(a) * f(mid) <= 0) b = mid; // root in left half else a = mid; // root in right half } return (a + b) / 2.0; } int main(void) { double root = bisection(1.0, 2.0, 1e-10); printf("root of x^3 - x - 1 = 0: %.10f\n", root); printf("check: f(root) = %e\n", f(root)); return 0; }
yn+1 = yn + h · f(xn, yn)k1 = f(xn, yn)k2 = f(xn+h, yn+h·k1)yn+1 = yn + h · (k1 + k2) / 2#include <stdio.h> #include <math.h> // dy/dx = f(x, y) double f(double x, double y) { (void)x; return -y; } double euler(double x0, double y0, double x_end, double h) { double x = x0, y = y0; while (x < x_end) { y = y + h * f(x, y); x = x + h; } return y; } double heun(double x0, double y0, double x_end, double h) { double x = x0, y = y0; while (x < x_end) { double k1 = f(x, y); double k2 = f(x + h, y + h * k1); y = y + h * (k1 + k2) / 2; x = x + h; } return y; } double rk4(double x0, double y0, double x_end, double h) { double x = x0, y = y0; while (x < x_end) { double k1 = f(x, y); double k2 = f(x + h/2, y + h/2 * k1); double k3 = f(x + h/2, y + h/2 * k2); double k4 = f(x + h, y + h * k3); y = y + h / 6 * (k1 + 2*k2 + 2*k3 + k4); x = x + h; } return y; } int main(void) { double exact = exp(-5.0); printf("exact y(5) = %.8f\n\n", exact); printf(" h | Euler | Heun | RK4\n"); printf("--------+-------------+-------------+-------------\n"); double hs[] = {0.5, 0.1, 0.01}; for (int i = 0; i < 3; i++) { double h = hs[i]; printf(" %5.3f | %.8f | %.8f | %.8f\n", h, euler(0, 1, 5, h), heun(0, 1, 5, h), rk4(0, 1, 5, h)); } return 0; }
m·x'' = -k·x by rewriting it as two first-order equations:dx/dt = v, dv/dt = -(k/m)·x#include <stdio.h> #define K 1.0 #define M 1.0 // right-hand side: returns dx/dt and dv/dt void deriv(double x, double v, double *dx, double *dv) { *dx = v; *dv = -(K / M) * x; } void rk4_step(double *x, double *v, double h) { double k1x, k1v, k2x, k2v, k3x, k3v, k4x, k4v; deriv(*x, *v, &k1x, &k1v); deriv(*x + h/2*k1x, *v + h/2*k1v, &k2x, &k2v); deriv(*x + h/2*k2x, *v + h/2*k2v, &k3x, &k3v); deriv(*x + h*k3x, *v + h*k3v, &k4x, &k4v); *x += h/6 * (k1x + 2*k2x + 2*k3x + k4x); *v += h/6 * (k1v + 2*k2v + 2*k3v + k4v); } int main(void) { double x = 1.0, v = 0.0; // amplitude 1, released from rest double h = 0.05; printf("t,x,v\n"); for (double t = 0; t <= 20; t += h) { printf("%.2f,%.4f,%.4f\n", t, x, v); rk4_step(&x, &v, h); } return 0; }
Check your understanding of this lesson!
double calculations?Values like 0.1 are repeating fractions in binary, so they're stored as approximations and introduce rounding error.
The x-intercept of the tangent line at the current estimate becomes the next estimate, rapidly converging toward the solution.
Kahan summation is a classic technique that tracks the low-order bits lost during floating-point addition and folds them back in as a compensation term.