Graham Gill

24 Reputation

2 Badges

18 years, 328 days

MaplePrimes Activity


These are answers submitted by Graham Gill

Thank you for your comments Robert and Acer.

Yes I do need the function to take on values other than the matrix entries M[i,j]. I need to be able to interpolate values. Currently I'm just using simple linear interpolation, but this might change. Yes, there is absolutely no need to have a table of ten separate functions f[1],...,f[10]. A single function f(n,x) with 1<=n<=10 an integer, x in [1,18] real, would do just as well. It's just because of math: I usually write f_n(x) instead of f(n,x) in such cases.

Acer's remark on how to try avoiding flattening of sequences inside seq() solved the whole problem for me. All I needed was the correct application of unevaluation quotes, and the remark showed me the way.

As I expected, I can do it all in one line:

f := (i,x) -> piecewise( seq( 'j - 1 < x <= j, (j - x) * `if`(j = 1, 0, M[i,j - 1]) + (x - j + 1) * M[i,j]', j = 1..18 ) );

[Trying to create a table g of functions g[1] to g[10] in a single statement this way seems quite a bit more complicated unfortunately. Again it boils down to getting just the right unevaluations I think - I'm not sure because I haven't got it working this way. The following does not work correctly:

g := 'g'; assign( seq( g[i] = proc(x) piecewise( seq( 'j - 1 < x <= j, (j - x) * `if`(j = 1, 0, M[i,j - 1]) + (x - j + 1) * M[i,j]', j = 1..18 ) ) end proc, i = 1..10 ) );

(The initial unassignment of g is necessary to avoid even more complicated unevaluations inside the assign(), in the case that g has already been assigned, such as when testing out variations on such a statement repeatedly!) This statement does not produce the desired effect because the i inside the unevaluation quotes in the inner seq() remains in the proc() as the name i with global scope.]

The test

seq( seq( evalb( f(i,j) = M[i,j]), j = 1..18), i = 1..10 );

returns a sequence of 180 trues, as desired, and statements such as

f(3,5.5);

do the right thing. However, the pretty-printed version of f produced when f is defined or produced by op(f) does not appear correct, or is at least confusing. op(f) produces a case defined function with a single case, "otherwise", the case definition being the literal contents of the call to piecewise() in the definition of f above.

Thanks for your comments. It's helpful to know that fsolve() uses Newton's method. Yes I guess that changing Digits would have the same effect as scaling up the unknowns (except possibly on execution time, but not in this problem once I'd broken it down - even changing Digits from 10 to 100 didn't make me wait much longer). However, going back to how I'd originally coded fsolve() when I got the floating point errors, in order to test the difference between changing Digits and scaling variables, is not something I want to revisit at the moment.

The data sets are floating point. 14 of the 18 equations are linear combinations of exponential functions with constant (floating point) coefficients, and the unknowns occur as arguments to the exponential functions. These 14 equations contain all 18 unknowns. The remaining 4 equations are simple linear interpolations to round out the system so that it will have a numerical solution. There is a certain amount of liberty in the problem in choosing the interpolations. Certain choices lead to an inconsistent system. The most desirable choice from my perspective turned out not to be desirable, as it lead to one unknown taking on a nonsensical value (in terms of the real world application), outside the range 0.02..0.05. Fortunately there was an alternative choice of interpolation which, although no more "correct" in the model than the other, at least led to sensible values for all unknowns. (The application is manually constructing yield curves for each day for a set of bonds over ten days. The bond data themselves were taken from real markets, so there was no "nice" solution to be expected.)

Yes I had to separate and reduce the system by hand as far as I could see, without more programming than the exercise is worth giving. (Already I've taken the exercise quite a bit farther than I needed to.) The choice of interpolating equations was necessary to do manually. The separating and reducing I did I was able to decide on manually without much work, but certainly solve() could help. Groebner wouldn't help here.

I could also have treated the entire system as linear with constant coefficients and used linear algebra to find the unknowns, then solved uniquely for the arguments to the exponential functions (provided the unknowns all turned out positive... which it's clear from the fsolve() solution they would have). However, I was seduced down the path of fsolve() and rather than abandon it wanted to know how I could help fsolve(), for future problems.

This kind of problem is more suited to MATLAB I suppose, but I only have the student version of Maple, not MATLAB, at the moment.

Page 1 of 1