#!/usr/bin/octave --silent unit_init(1, {}); unit_test_equals("sin[0] == 0", ... 0, ... divided_difference(@sin, 0)); unit_test_equals("sin[0, pi] == 0", ... 0, ... divided_difference(@sin, [0,pi])); unit_test_equals("sin[0, pi, 2*pi] == 0", ... 0, ... divided_difference(@sin, [0,pi,2*pi])); unit_test_equals("zero order divided_difference_coefficients", ... [1], ... divided_difference_coefficients([0])); unit_test_equals("first order divided_difference_coefficients", ... [-1, 1] / pi, ... divided_difference_coefficients([0, pi])); unit_test_equals("second order divided_difference_coefficients", ... [1, -2, 1] / (2*pi^2), ... divided_difference_coefficients([0, pi, 2*pi])); unit_test_equals("1 is odd", ... true, ... odd(1)); unit_test_equals("1 is not even", ... false, ... even(1)); unit_test_equals("2 is not odd", ... false, ... odd(2)); unit_test_equals("2 is even", ... true, ... even(2)); expected_A = [1, 0, 0, 0, 0; ... 16, -32, 16, 0, 0; ... 0, 16, -32, 16, 0; ... 0, 0, 16, -32, 16; ... 0, 0, 0, 0, 1]; unit_test_equals("Homework #1 problem #1 Poisson matrix is correct", ... true, ... expected_A == poisson_matrix(4, 0, 1)); g = @(x) 1 + atan(x); expected_fp = 2.1323; tol = 1 / 10^10; x0 = 2.4; unit_test_equals("Homework #2 problem #5 fixed point is correct", ... expected_fp, ... fixed_point_method(g, tol, x0)); h = 0.5; g1 = @(u) 1 + h*exp(-u(1)^2)/(1+u(2)^2); g2 = @(u) 0.5 + h*atan(u(1)^2 + u(2)^2); my_g = @(u) [g1(u), g2(u)]; tol = 1 / 10^9; u0 = [1,1]; expected_fp = [1.0729, 1.0821]; unit_test_equals("Homework #3 problem #3i fixed point is correct", ... expected_fp, ... fixed_point_method(my_g, tol, u0)); f = @(x) x^6 - x - 1; f_prime = @(x) 6*x^5 - 1; tol = 1/1000000; x0 = 2; expected_root = 1.1347; unit_test_equals("Newton's method agrees with Haskell", ... expected_root, ... newtons_method(f, f_prime, tol, x0)); f1 = @(u) u(1)^2 + u(1)*u(2)^3 - 9; f2 = @(u) 3*u(1)^2*u(2) - u(2)^3 - 4; f = @(u) [f1(u); f2(u)]; ## The partials for the Jacobian. f1x = @(u) 2*u(1) + u(2)^3; f1y = @(u) 3*u(1)*u(2)^2; f2x = @(u) 6*u(1)*u(2); f2y = @(u) 3*u(1)^2 - 3*u(2)^2; ## f_prime == Jacobian. f_prime = @(u) [ f1x(u), f1y(u); f2x(u), f2y(u) ]; tol = 1 / 10^12; u0 = [1.2; 2.5]; expected_root = [1.33635; 1.75424]; [actual_root, iterations] = newtons_method(f, f_prime, tol, u0); unit_test_equals("Homework #3 problem #4 root is correct", ... expected_root, ... actual_root);