Combining Central Difference Scheme and Gaussian Elimination to Solve Matrix

One part of a question requires me to discretize the following equation, using the central difference method. -u''=x-(1/2), for x ∈ [0, 1], u(0) =2, u(1)=-1

I've never taken Numerical methods, so I'm having trouble turning this into a matrix problem basically. I have a Gaussian Elimination routine ready to solve for x, but I'm still confused how to apply the Numerical methods techniques correctly to create "u" and "b" matrices for the equation AU=b.

Right now I have the following code, and perhaps someone can send me in the right direction from here.

%Boundary Conditions
x_0=0;
x_n=1;
b_0=2;
b_n=-1;

%Other variables
dx = 0.1;
x=[x_0:dx:x_n]; %Initializing x vector
n = length(x); %Number of discrete points calculated
h = dx; %Distance between points
u = zeros(1,n); 
A = zeros(n,n);

%Adding Bound Conditions to arrays
x(1)=x_0;
x(n)=x_n;
b(1)=b_0;
b(n)=b_n;

%Create "A" Matrix
for i=1:n-2,
A(i,i) = -2/h^2; 
A(i+1,i) = 1/h^2; 
A(i,i+1)= 1/h^2;
end
A(n-1,n-1) = -2/h^2;

f = inline('x-0.5','x'); 
f(x(i))=(u(i-1)-2*u(i)+u(i+1))/h^2;

for j = 2:n-1 
x(j) = x(1)+(j-1)*h;
F(j) = (u(i-1)-2*u(i)+u(i+1))/h^2;
end;

Answer 1

[Zoltán Csáti]

Your Gauss-elimination program takes effect after this line:

%%Solve the linear system

If you must use your Gauss-elimination method instead of the built-in method, use your code instead of u = A\b. If it works for the 3x3 magic example, it will probably work for this task too since the condition number of the coefficient matrix A is low.

To answer your second question: Since the endpoint values are known from the boundary conditions, it is enough to solve for the inner values. Later it is augmented with the known values, u_0 and u_n.

[SB]

Edit: I figured out what was wrong with the code. The first issue was re-using "x" in the Gaussian Elimination part, which caused issues when trying to plot against the original "x" at the end. I was on the right track with the use of the variables in GE, but made some tweaks and finally got the code working properly. Final result is below.

- - - - - - - - - - - - - - - - - - - - - - - - - - -

Ok, so my specific question for that final part is how to assign the variables properly to perform Gaussian Elimination. My plot results didn't look right to me.

The combined code basically assigns U=A (for LU factorization), does not calculate "u" from central difference scheme, but does calculate "x" as part of the LU factorization. So, my final result is "x", not "u". Am I missing something here? I use the same "b" for both parts. Code is below.

Also, where does the formula for the central difference approximation come in? In my notes, it's written as follows, for the 2nd derivative: u(x)=(u(x+h)-2*u(x)+u(x-h))/h^2

Thank you again for your time!!

      %%Centered Difference Scheme
    % Boundary Conditions
    x_0 = 0;
    x_n = 1;
    u_0 = 2;
    u_n = -1;

    n = 20; % Number of discrete points calculated
    h = (x_n-x_0)/(n-1); % Distance between points
    x = (x_0:h:x_n)'; % Initializing x vector

    % Other variables

    % Create the tri-diagonal coefficient matrix "A"
    A = 1/h^2*(diag(-2*ones(1,n-2)) + diag(ones(1,n-3),1) + diag(ones(1,n-3),-1));

    % Create "b" matrix
    f = 1/2-x(2:n-1);
    b = f; 
    b(1) = f(1)-u_0/(h^2); 
    b(n-2) = f(n-2)-u_n/(h^2);

    %Exact solution
    hold on
    u_exact = @(x2) (1/4)*x2.^2-(1/6)*x2.^3+(-73/24)*x2+2;
    fplot(u_exact,[-1 1 -1 2],'LineWidth',2,'r');
    %%Gaussian Elimination 
    % Naive Version without Pivoting
    %  A=magic(3);  % Values for testing algorithm
    %  b=[1;3;5];

    [m,m]= size(A); %Assuming square matrix here
    U=A;
    l = zeros(m,m-1);
    b2=b;
    % Working Forward
    for k = 1:m-1 
        for j = k+1:m 
            l(j,k) = U(j,k)/U(k,k); %Divide by element above
            U(j,k:m) = U(j,k:m)-l(j,k)*U(k,k:m); %Subtracting row above, multipled by L
            b2(j) = b2(j) - l(j,k)*b2(k); %Applying the same math to the b side
        end
    end

    % Working Backward to Solve for X
    x_GE = zeros(m,1);
    x_GE(m) = b2(m)/U(m,m);

    for k = m-1:-1:1
        x_GE(k) = (b2(k) - U(k,k+1:m)*x_GE(k+1:m))/U(k,k);
    end
    x_GE = [u_0; x_GE; u_n];

    % %Plot the No-Pivot Solution
    % figure;
    % plot(x,b2,'o-')

    %Numerical Solution

    hold on
    plot(x,x_GE,'x-','LineWidth',0.001)
    %

[Zoltán Csáti]

I recommend you to separate the FD approximation from the linear system solving. Your Gauss elimination code must work in all cases. However it does not work for me. So the "Gaussian Elimination" block should contain the solver call: u = gausselim(A,b); where gausselim is your function which solves the discretized system.

"Also, where does the formula for the central difference approximation come in? In my notes, it's written as follows, for the 2nd derivative: u(x)=(u(x+h)-2*u(x)+u(x-h))/h^2". It is the same formula that I used. Note that u(x_i) = u_i, u(x_i+h) = u_(i+1) and u(x_i-h) = u_(i-1), i=2,...,n-1. Now you form the matrix vector product A*u, where u = (u_2;u_3;...;u_(n-1)) and A is a tridiagonal matrix (because for each i, there is u_(i-1), u_i and u_(i+1)).

Answer 2

[Zoltán Csáti]

I preserved the structure of your code, but modified it. Now it perfectly works.

%%Boundary Conditions
x_0 = 0;
x_n = 1;
u_0 = 2;
u_n = -1;

%%Other variables
h = 0.2; % Distance between points
x = (x_0:h:x_n)'; % Initializing x vector
n = length(x); % Number of discrete points calculated
% x = x(2:n-1); % Extract the inner points

%%Create the tridiagonal coefficient matrix
A = 1/h^2*(diag(-2*ones(1,n-2)) + diag(ones(1,n-3),1) + diag(ones(1,n-3),-1));

%%Create the right hand side
f = 1/2 - x(2:n-1);
b = f; b(1) = f(1)-u_0/(h^2); b(n-2) = f(n-2)-u_n/(h^2);

%%Solve the linear system
u = A\b;
u = [u_0; u; u_n];

%%Plot the discrete solution
figure;
plot(x,u,'o-')

%%Compare it with the exact solution
hold on;
u_e = @(x) 95/48-(x-1/2).^3/6-(71*x)/24;
fplot(u_e,[0 1 -1 2],'r');

%%Error at the nodes
err = u_e(x) - u;

[SB]

Thank you! I suppose a part of my question still remains though: how do I tie in my Gaussian Elimination code to this? Basically I'm supposed to use that central difference technique with Gaussian elimination.

My GE code with no pivoting is below, as an example.

%%Gaussian Elimination - Without Pivoting
 %A=magic(3);  % Values for testing algorithm
 %b=[1;3;5];
u=A;
[m,m]= size(A); %Assuming square matrix here
l = zeros(m,m-1);
% Working Forward
for k = 1:m-1 
    for j = k+1:m 
        l(j,k) = u(j,k)/u(k,k); %Divide by element above
        u(j,k:m) = u(j,k:m)-l(j,k)*u(k,k:m); %Subtracting row above, multipled by L
        b(j) = b(j) - l(j,k)*b(k); %Applying the same math to the b side
    end
end

% Working Backward to Solve for X
x = zeros(m,1);
x(m) = b(m)/u(m,m);

for k = m-1:-1:1
    x(k) = (b(k) - u(k,k+1:m)*x(k+1:m))/u(k,k);
end

[SB]

Also, should "A" be the same dimensions as "n"? Your code for A gives an n-2 by n-2 matrix.