Combining Central Difference Scheme and Gaussian Elimination to Solve Matrix
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
SB
am 19 Okt. 2014
Kommentiert: Zoltán Csáti
am 22 Okt. 2014
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 Kommentare
Akzeptierte Antwort
Zoltán Csáti
am 21 Okt. 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 Kommentare
Zoltán Csáti
am 22 Okt. 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)).
Weitere Antworten (1)
Zoltán Csáti
am 20 Okt. 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;
Siehe auch
Kategorien
Mehr zu Operating on Diagonal Matrices finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!