#!/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);



f = @(x,y) y;
x0 = 0;
y0 = 1;
h = 1;
actual_y = forward_euler1(x0, y0, f, h);
expected_y = 2;

unit_test_equals("Forward Euler works for one step", ...
		 expected_y, ...
		 actual_y);
