Problem when solving the Thomas algorithm

[JohanCarlsen]

Hi, I'm stuck. I can't seem to find where the error is in my code.

`// initial values double x_0 = step; double u_0 = f(x_0); double b_0 = b; double g_0 = u_0;

// initialize vectors x(0) = x_0; u(0) = u_0; b_tilde(0) = b_0; g_tilde(0) = g_0;

// forward substitution for (int i = 0; i < n; i++) { x(i+1) = x(i) + step; u(i+1) = f(x(i+1)); b_tilde(i+1) = b - a / b_tilde(i) c; g_tilde(i+1) = u(i+1) - a / b_tilde(i) g_tilde(i); }

v(n) = g_tilde(n) / b_tilde(n);

// backwards substitution for (int i = n-1; i >= 0; i--) { v(i) = (g_tilde(i) - c * v(i+1)) / b_tilde(i); }`

This produces as plot, which for u(x) seems ok, but for v(x) is very strange. Does anyone spot the error, or am I forever doomed?

[fridasor]

I may be mistaken, but it looks like maybe you've forgotten to multiply f(x+1) by h^2?

[JohanCarlsen]

Oh! But aren’t g_i just the discretized version of f(x)? If not, I may have to rethink a lot..

[anderkve]

Hi @JohanCarlsen!

I only have a couple of minutes just now, so this reply will be very quick. A few things:

• If you have a version of this code on GitHub, it may be easier to provide a link to the code directly, since then we can see how everything is exactly defined, e.g. the variables a, b, c, etc.

• No, g_i is not just the discretized version of f(x). As @fridasor mentions, there is a factor of h^2 involved. (And in the general case, with non-zero boundary conditions, there is also additional terms for the first and last g element.) I'd recommend going back to Problem 4 and make sure you see how this h^2 factor appears, before working more on the code.

• On the code side: It looks like you are trying to do many things at once in your code, i.e. both solving the matrix equation A v = g with Gaussian elimination, and filling vectors x and u at the same time. But the vectors x and u should not be involved at all in the Gaussian elimination algorithm for solving A v = g. So I'd recommend dividing up your code so that one part of the code (e.g. a separate function) takes care of only the problem of solving A v = g. (For instance, I see that u(i+1) appears in your forward substitution part, which it certainly shouldn't, so splitting apart your code will probably help you avoid mistakes like that.)

Good luck with the debugging! :)

[mikkelme]

Hi @JohanCarlsen,

First of all, I agree with the advice given already. The general mistakes seem to be that you confuse u and g and also have some inconsistencies with the indexing. As @anderkve mentioned, it is a bit difficult to debug without seeing how you defined the variables and the length of your arrays, but I made an attempt regardless. Bear in mind that it really makes things more complicated when you fill in arrays and solve the matrix in the same loops (as mentioned already), but for educational purposes I kept that in the debugged version for better transparency of what needs to be changed. Also, in your final version a and c should be in array form (and not just a single variable) such that you have more flexibility for different tridiagonal matrices, but I did not change that here.

My assumption is that your vector v for the inner solution has length n. The following changes should solve the problems.

// Change notation 
double h = step;
double hh = h*h; 

// Initialize vectors
x(0) = h;
b_tilde(0]) = b;
g_tilde(0) = f( x(0) )*hh;

// Forward substitution
for (int i = 1; i < n; i++) // Changed starting point to 1 and indexing from, i+1 to i, inside the loop
{
    // Fill in arrays (should be done in a separate loop for better readability)
    x(i) = x(i-1) + h; 
    g(i) = f( x(i) ) * hh;

    // Actual forward substitution part 
    b_tilde(i) = b - a / b_tilde[i-1] * c;
    g_tilde(i) = g(i) - a / b_tilde(i-1) * g_tilde(i-1);
}

// Backwards substitution
v(n-1) = g_tilde(n-1) / b_tilde(n-1);   // changed index from n to n-1
for (int i = n-2; i >= 0; i--)   // changed start point to n-2
{
    v(i) = (g_tilde(i) - c * v(i+1)) / b_tilde(i);
}

// End of code

Now the inner part of the solution should be contained in x, v such that the full solution is given as x = [0, x(0), x(1), ..., x(n-1), 1] v = [0, v(0), v(1), ..., v(n-1), 0]

[JohanCarlsen]

Thank you all. It worked!