%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.

Google Online Preview   Download