Varijacione metode (Variational Methods) by Mihailo Despotovic, September 1995

My seminar paper for the Numerical Methods class of professor Desa Radunović, during my undergraduate studies of Mathematics and Computer Science at the University of Belgrade, Faculty of Mathematics. It solves the boundary value problem p(x)u″ + q(x)u′ + r(x)u = f(x), u(a) = c, u(b) = d, approximately, by a polynomial of degree n that meets the two boundary conditions, found by three variational methods: Galerkin's, collocation, and least squares. It draws the result in VGA graphics and, given the exact solution, draws that too and measures the error for every degree from 2 to 8. The functions are typed in as text, like sin(x)+x*x. It runs here as it did in DOS, in Serbian, with the English translation beside it. Click the screen and type, or use the buttons.

The menu (MENU.C) and the expression evaluator (EVALL.C) are Vuksan Pejović's libraries, as is the symbolic differentiation in DIFF/, which an earlier version of the program used.

The disk

P saves the equation to a file and loads it back; after a solution, D saves it to a file of results, adding to it. The files are kept here while the page is open.

How it works

The solution is sought as u = w + Σ c_kφ_k, k = 2 … n. Here w = −koef1 − koef2·x is the straight line through the two boundary values, and φ_k(x) = x^k−2(x − a)(x − b) are zero at both ends, so any c_k keep the boundary conditions. Putting u into the equation leaves a residual, and the three methods choose the c_k differently: Galerkin's makes the residual orthogonal to every φ_j; collocation makes it zero at n − 1 points spread over [a, b]; least squares makes the integral of its square as small as possible. Each gives n − 1 linear equations. The program builds the functions it integrates as text, as ((1)*2+(2)*(2*x-1)+(3)*(x*x-x))*(x*x-x), has Vuksan's evaluator parse them, integrates them by Simpson's rule, doubling the pieces until the estimated error is below 0.00001, and solves the equations by Gaussian elimination.

The bugs, and the fixes

• The equations had one term too few. Putting u = w + Σ c_kφ_k into the equation moves L(w) = q·w′ + r·w to the right, and r·w = −r·koef1 − r·koef2·x; the program left out the last part. So whenever r(x) isn't 0 and the two boundary values differ, all three methods solved a different equation. For the default one, u″ + 2u′ + 3u = 0, u(0) = 0, u(1) = 1, they all settled on about 2.78x − 2.78x² + …, where the true solution starts 3.89x − 3.89x² + …; at x = 0.5, Galerkin with n = 2 gave 0.86, and the exact solution is 1.08. The small parameter's equation has r = 0, so it came out right.
• Vuksan's evaluator read a − b − c as a − (b − c), and a / b / c as a / (b / c). The program writes its basis functions with that pattern on some intervals: on [−1, 2], φ₂ is (x*x-x-2), read as x² − x + 2. Then every method is wrong: for u″ + u = 0, u(−1) = u(2) = 1, all three give about 1.6 at x = 0.5, where the exact solution is 14.14. The evaluator also can't take more than 31 numbers in one expression, since it stands each for a character below 32, fails a power of a power after anything longer than one character, as in (x+1)^2^3, and writes one byte past the memory it asks for.
• Loading the data didn't take the new boundary values in: the line w was computed when the equation was typed in, and not again after a load, so until you changed the equation, the old boundary values were used.
• The graph's height was measured on the wrong function. With an exact solution, its values were meant to widen the vertical range, but the program evaluated whatever expression had been parsed last, the last function integrated. So the exact solution could run off the top or the bottom of the screen.
• The error chart draws bars that shouldn't be there. A missing pair of braces in if (odstoj[n2-2]>0) makes only the first line conditional, so where the error is 0, or the system couldn't be solved, a bar is drawn anyway, with the height of the one before; the first one with a height never set at all. The errors are also capped at 10, to keep the chart readable, for every n but 2.
• The largest error could miss the end of the interval. It is sought by adding 1/633 of [a, b] at a time, and whether the sum reaches b exactly depends on the rounding. With Borland C's 80-bit numbers it did, for [0, 1].
• Smaller ones: the “exact” solution for the small parameter uses 2.718 for e; loading a 39-character expression reads it in two pieces, which shifts every line after it, and the numbers are saved with 6 decimals; and the program also wrote u(x) in terms of the basis, through a pointer that points nowhere, and never showed it.

Fixed corrects all of these, with a − b − c and a / b / c read left to right, any number of constants, and exp(x) in the small parameter's solution. At the higher degrees, Galerkin's method can still go astray, as the basis functions grow nearly alike; that is the method with this basis, not a bug.

The port was checked against the original C, compiled with a simulated DOS screen, keyboard and BGI, in 112 sessions: every method and degree, equations, intervals, exact solutions, the error charts, saving and loading. In 109 the screens, the drawing and the files are the same. Borland C computed in 80 bits and JavaScript computes in 64, so the C was also compiled with 80-bit numbers. The one place where that changed what the program shows is the largest error above, and the page does that in 80 bits too; that accounts for two of the other three sessions. In the third, five pixels of a curve fall exactly between two rows. The 80-bit build differs besides only in the bars drawn with a height that was never set.