Combining Central Difference Scheme and Gaussian Elimination to Solve Matrix
2 views (last 30 days)
Show older comments
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;
0 Comments
Accepted Answer
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
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)).
More Answers (1)
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;
See Also
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!