Combining Central Difference Scheme and Gaussian Elimination to Solve Matrix

9 views (last 30 days)
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;

Accepted Answer

Zoltán Csáti
Zoltán Csáti on 21 Oct 2014
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.
  2 Comments
SB
SB on 21 Oct 2014
Edited: SB on 22 Oct 2014
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
Zoltán Csáti on 22 Oct 2014
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)).

Sign in to comment.

More Answers (1)

Zoltán Csáti
Zoltán Csáti on 20 Oct 2014
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;
  2 Comments
SB
SB on 20 Oct 2014
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
SB on 20 Oct 2014
Also, should "A" be the same dimensions as "n"? Your code for A gives an n-2 by n-2 matrix.

Sign in to comment.

Categories

Find more on Operating on Diagonal Matrices in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!