Add the is_upper_triangular() function and its tests.
[octave.git] / forward_euler.m
1 function coefficients = forward_euler(integer_order, xs, x)
2 ##
3 ## Return the coefficients of u(x0), u(x1), ..., u(xn) as a vector.
4 ## Take for example a first order approximation, with,
5 ##
6 ## xs = [x0,x1,x2,x3,x4]
7 ##
8 ## f'(x1) ~= [f(x2)-f(x1)]/(x2-x1)
9 ##
10 ## This would return [0, -1/(x2-x1), 2/(x2-x1), 0, 0]. This aids the
11 ## solution of linear systems.
12 ##
13 ##
14 ## INPUTS:
15 ##
16 ## * ``integer_order`` - The order of the derivative which we're
17 ## approximating.
18 ##
19 ## * ``xs`` - The vector of x-coordinates.
20 ##
21 ## * ``x`` - The point `x` at which you'd like to evaluate the
22 ## derivative of the specified `integer_order`. This should be an
23 ## element of `xs`.
24 ##
25 ##
26 ## OUTPUTS:
27 ##
28 ## * ``coefficients`` - The vector of coefficients, in order, of
29 ## f(x0), f(x1), ..., f(xn).
30 ##
31
32 if (integer_order < 0)
33 ## You have made a grave mistake.
34 coefficients = NA;
35 return;
36 end
37
38 if (integer_order == 0)
39 coefficients = x;
40 return;
41 end
42
43 if (length(xs) < 2)
44 ## You can't approximate a derivative of order greater than zero
45 ## with zero or one points!
46 coefficients = NA
47 return;
48 end
49
50 if (even(integer_order))
51 offset_b = integer_order / 2;
52 offset_f = offset_b;
53 else
54 ## When the order is odd, we need one more "forward" point than we
55 ## do "backward" points.
56 offset_b = (integer_order - 1) / 2;
57 offset_f = offset_b + 1;
58 end
59
60 ## Zero out the coefficients for terms that won't appear. We compute
61 ## where `x` is, and we just computed how far back/forward we need to
62 ## look from `x`, so we just need to make the rest zeros.
63 x_idx = find(xs == x);
64 first_nonzero_idx = x_idx - offset_b;
65 last_nonzero_idx = x_idx + offset_f;
66 leading_zero_count = first_nonzero_idx - 1;
67 leading_zeros = zeros(1, leading_zero_count);
68 trailing_zero_count = length(xs) - last_nonzero_idx;
69 trailing_zeros = zeros(1, trailing_zero_count);
70
71 targets = xs(first_nonzero_idx : last_nonzero_idx);
72
73 # The multiplier comes from the Taylor expansion.
74 multiplier = factorial(integer_order);
75 cs = divided_difference_coefficients(targets) * multiplier;
76
77 coefficients = horzcat(leading_zeros, cs, trailing_zeros);
78 end