%Direct Solve by Gaussian Elimination
BVP Sample MATLAB Codes
1. %Direct Solve by Gaussian Elimination
clear; close all; clc; format short
for dx = [1/2 1/3 1/5 1/100]
k = 0;
le = 0; %left endpoint
re = 1; %right endpoint
n = round((re-le)/dx)+1;
x = linspace(le,re,n);
a = zeros(1,n); b=a;d=a;c=a; residual = a; y=a;
A = zeros(n);
i = 1;
a(i) = 0;
b(i) = 1;
c(i) = 0;
d(i) = 0;
for i = 2:n-1
a(i) = 1/dx^2;
b(i) = -2/dx^2;
c(i) = 1/dx^2;
d(i) = 1;
end
i = n;
a(i) = 0;
b(i) = 1;
c(i) = 0;
d(i) = 0;
%SOLVE TRI-DIAGONAL LINEAR SYSTEM
bb = zeros(1,n);
db = bb;
gs = db;
bb(1)=b(1);
db(1)=d(1);
for j=2:n;
bb(j)=b(j)-(a(j)/bb(j-1))*c(j-1);
db(j)=d(j)-(a(j)/bb(j-1))*db(j-1);
end
gs(n)=db(n)/bb(n);
for j=n-1:-1:1;
gs(j)=(db(j)-c(j)*gs(j+1))/bb(j);
end
y = gs;
for i = 2:n-1
residual(i) = 1/dx^2*(y(i+1)-2*y(i)+y(i-1))-1;
end
error = max(abs(residual));
plot(le:dx:re,y,'k'),hold on%,
end
title('Beam Bending Problem, u"=1, u(0)=0, u(1)=0')
2.
%Iterative Methods: 'Jacobi','Gauss-Seidel','SOR'
%
%We will solve the BVP:
% y''+y = 1, y(-pi/4) = 0, y(pi/4) = 0
%With many step sizes
clear %clear all variables in memory
close all %close all open windows (plots, etc...)
clc %clear the Command Window (just nice to have a clean slate)
format short e %can declare a format for output if you want to
TOL = 1e-6; %set tolerance
le = 0; %left endpoint
re = 1; %right endpoint
ik = 0;
disp('TOL = 1e-6')
disp('N JACOBI GAUSS-SEIDEL SOR')
for dx = [1/2 1/3 1/5 1/100], ik = ik+1;
n = (re-le)/dx+1; %size of problem
x = linspace(le,re,n); %fast way to make the x-vector (for plotting)
scheme={'Jacobi','Gauss-Seidel','SOR'}; %use this for plotting titles
omega = 2/sqrt(1+sin(pi*dx)); %set optimal omega
for method = 1:3
%Re-Initialize all parameters for each method
y = zeros(1,n); ynew = y; residual = y;
k = 0; error = 1;
while error > TOL && k < 1e6
if method == 1 %Jacobi
for i = 2:n-1 % - lags the updated values to the next iteration
ynew(i) = 1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1);
end
y = ynew;
elseif method == 2 %Gauss-Seidel
for i = 2:n-1 % - uses the updated values as soon as available
y(i) = 1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1);
end
elseif method == 3 %SOR
for i = 2:n-1 % - hybrid Gauss-Seidel: x(i+1) = x(i)+omega*res(i)
y(i) = (1-omega)*y(i) ...
+ omega *(1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1));
end
end
%Check the residual to see how close we are to a solution
for i = 2:n-1
residual(i) = 1/dx^2*(y(i+1)-2*y(i)+y(i-1))-1;
end
error = max(abs(residual)); %set error
k = k+1; %update count
end %while
count(method) = k;
%figure(ik)
subplot(3,1,method),
plot(x,y), hold on, ylim([-.15 0])
title(scheme(method));
end
fprintf('%5.0f %5.0f %5.0f %5.0f\n',n, count(1),count(2),count(3))
end
3. Show difference between residual and iterations between the methods on a single figure
%Iterative Methods: 'Jacobi','Gauss-Seidel','SOR'
%
%We will solve the BVP:
% y''+y = 1, y(-pi/4) = 0, y(pi/4) = 0
clear %clear all variables in memory
close all %close all open windows (plots, etc...)
clc %clear the Command Window (just nice to have a clean slate)
format short e %can declare a format for output if you want to
TOL = 1e-6; %set tolerance
le = 0; %left endpoint
re = 1; %right endpoint
ik = 0;
disp('TOL = 1e-6')
disp('N JACOBI GAUSS-SEIDEL SOR')
error = zeros(700,3); error(1,1:3)=1;
for dx =1/10,
n = (re-le)/dx+1; %size of problem
x = linspace(le,re,n); %fast way to make the x-vector (for plotting)
scheme={'Jacobi','Gauss-Seidel','SOR'}; %use this for plotting titles
omega = 2/(1+sin(pi*dx)); %set optimal omega
for method = 1:3,ik= ik+1;
%Re-Initialize all parameters for each method
y = zeros(1,n); ynew = y; residual = y;
k = 1; %error = 1;
while error(k,ik) > TOL && k < 1e6
if method == 1 %Jacobi
for i = 2:n-1 % - lags the updated values to the next iteration
ynew(i) = 1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1);
end
y = ynew;
elseif method == 2 %Gauss-Seidel
for i = 2:n-1 % - uses the updated values as soon as available
y(i) = 1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1);
end
elseif method == 3 %SOR
for i = 2:n-1 % - hybrid Gauss-Seidel: x(i+1) = x(i)+omega*res(i)
y(i) = (1-omega)*y(i) ...
+ omega *(1/(-2/dx^2)*(-1/dx^2*(y(i+1)+y(i-1))+1));
end
end
%Check the residual to see how close we are to a solution
for i = 2:n-1
residual(i) = 1/dx^2*(y(i+1)-2*y(i)+y(i-1))-1;
end
k = k+1; %update count
error(k,ik) = max(abs(residual)); %set error
end %while
count(method) = k;
end
fprintf('%5.0f %5.0f %5.0f %5.0f\n',n, count(1),count(2),count(3))
end
semilogy(1:700,error(:,1:3)),title('dx = [1/10]'),xlabel('Iterations'),ylabel('Residual')
legend('Jacobi','Gauss-Seidel','SOR')
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- solve by elimination calculator
- gaussian elimination matrix
- gaussian elimination rules
- gaussian elimination calculator
- gaussian elimination method
- gaussian elimination 2x2 matrix
- the gaussian elimination method
- solve matrix by gaussian elimination
- gaussian elimination 3x3
- gaussian elimination explained
- gaussian elimination method example
- gaussian elimination method 4x4