# Numerical Analysis Lab

Document Sample

```					Numerical Analysis Lab

Lecture 4: Solutions of Equations in One Variable

Hongyu Liu

Autumn 2010
The Bisection Method

• To find a solution to f(x)=0 given the
continuous function f on the intervl [a,b],
where f(a) and f(b) have opposite signs.
Pseudo-code

Input endpoints a,b; TOL; maximum no. of iteration N0.
Output approx. sol. p or message of failure.
Step 1 Set i=1;
FA=f(a);
Step 2 While i<=N0 do Step 3—6.
Step 3 Set p=a+(b-a)/2;
FP=f(p).
Step 4 If FP=0 or (b-a)/2<TOL then
OUTPUT(p);
Stop.
Step 5 Set i=i+1.
Step 6 If FA*FP>0, then set a=p;
FA=FP
else set b=p.
Step 7 Output (‘Method failed after N0 iterations, N0=’, N0)
STOP
function[b]=f(x);
b=x^2-4;
a=-1;
b=4;
TOL=1e-04;
N0=30;
k=1;
FA=f(a);
while k<=N0;
p=a+(b-a)/2;
FP=f(p);
if abs(FP)<TOL | (b-a)/2<TOL;
p
break;
end
k=k+1;
if FA*FP>0;
a=p;
FA=FP;
else
b=p;
end
end
if k>30;
disp('Method Failed after N0 iterations');
N0
end
Excercise

1. Formulate the bisection method as a function:

Bisec(a,b,TOL,N0)
2. Use your program to determine a root for
f(x)=x^3+4x^2-10=0 in [1,2] that is accurate with
10^{-4}.
Fixed-Point Iteration

For a function g(x), find the point g(p)=p
Pseudo-code

Input initial approximation p0; tolerance TOL; maximum no. of iteration N0
Output approximate solution p or message of failure.
Step 1 Set k=1.
Step 2 While k<=N0 do Steps 3—6
Step 3 Set p=g(p0). (compute p_k)
Step 4 If |p-p0|<TOL then
OUTPUT p; (the procedure is successful)
STOP.
Step 5 Set k=k+1;
Step 6 Set p0=p. (Update p0)
Step 7 OUTPUT Message of Failure
STOP
p0=0;
TOL=1e-04;
N0=20;
k=1;
while k<=N0;
p=g(p0);
if abs(p-p0)<TOL;
p
break;
end
k=k+1;
p0=p;
end
if k>N0;
disp('The method failed to work after N0 iterations');
N0
end
Exercise

The equation x^3+4x^2-10 has a unique root in [1,2]. Reformulate the
equation as a fixed-point problem in the following five ways.
(a) x=x-x^3-4x^2+10
(b) x=(10/x-4x)^{1/2}
(c) X=1/2(10-x^3)^{1/2}
(d) X=(10/(4+x))^{1/2}
(e) X=x-(x^3+4x^2-10)/(3x^2+8x)

With an initial approximation p0=1.5, find the solution.
Newton’s Method

Newton’s method is one of the most powerful and
well-known numerical method for solving a root
-finding problem.
Pseudo-code

Input initial approximation p0; tolerance TOL; maximum no. of iteration N0
Output approximate solution p or message of failure.
Step 1 Set k=1.
Step 2 While k<=N0 do Steps 3—6
Step 3 Set p=p0-f(p0)/f’(p0). (compute p_k)
Step 4 If |p-p0|<TOL then
OUTPUT p; (the procedure is successful)
STOP.
Step 5 Set k=k+1;
Step 6 Set p0=p. (Update p0)
Step 7 OUTPUT Message of Failure
STOP
Exercise

Use Newton’s method to find solutions accurate to within 1e-05 to the
following problems
1.     1-4xcos(x)+2x^2+cos(2x), for x in [0,1]
2.    x^2-2xexp(-x)+exp(-2x)=0, for x in [0,1]
Secant Method
Pseudo-code

Input initial approximation p0, p1; tolerance TOL;
maximum no. of iteration N0
Output approximate solution p or message of failure.
Step 1 Set k=2; q0=f(p0); q1=f(p1).
Step 2 While k<=N0 do Steps 3—6
Step 3 Set p=p1-q1(p1-p0)/(q1-q0). (compute p_k)
Step 4 If |p-p1|<TOL then
OUTPUT p; (the procedure is successful)
STOP.
Step 5 Set k=k+1;
Step 6 Set p0=p1;q0=q1;p1=p;q1=f(p). (Update p0,q0,p1,q1)
Step 7 OUTPUT Message of Failure
STOP
Gaussian Elimination with Partial Pivoting

To solve the nxn linear system
Pseudo-code

Input number of unknowns and equations n; agumented matrix A=[aij] where 1<=i<=n
and 1<=j<=n+1.
Output solutions x1,…,xn or message that the linear system has no unique solution.
Step 1 For i=1,….,n set NROW(i)=i. (Initialize row pointer)
Step 2 For i=1,…,n-1 do Steps 3—6 (elimination process)
Step 3 Let p be the smallest integer with i<=p<=n and
|a(NROW(p),i)|=maxi<=j<=n |a(NROW(j),i)|.
Step 4 If a(NROW(p),i)=0 then OUTPUT(‘no unique solution’); STOP.
Step 5 If NROW(i)~=NROW(p) then set NCOPY=NROW(i);
NROW(i)=NROW(p);
NROW(p)=NCOPY.
(simulated row interchange)
Step 6 For j=i+1,…,n do Steps 7, 8.
Step 7 Set m(NROW(i),j)=a(NROW(j),i)/a(NROW(i),i).
Step 8 Perform (ENROW(j)-m(NROW(j),i)*ENROW(i))- ENROW(j)
Step 9 If a(NROW(n),n)=0, then OUTPUT(‘no unique solution’); STOP.
Step 10 Set xn=a(NROW(n),n+1)/a(NROW(n),n).
Step 11 For i=n-1,…,1

set

Step 12 OUTPUT x.
STOP
n=3;
A=[1 2 3 2;0 1 0 3;0 0 1 4];
for k=1:3;
NROW(k)=k;
end
for k=1:n-1;
[U,d]=max(A(k:n,k));
p=d+k-1;
if A(p,k)==0;
disp('no solution');
break;
end
if NROW(k)~=NROW(p);
NCOPY=NROW(k);
NROW(k)=NROW(p);
NROW(p)=NCOPY;
end
for l=k+1:n;
m=A(NROW(l),k)/A(NROW(k),k);
T=A(NROW(l),:)-m*A(NROW(k),:);
A(NROW(l),:)=T;
end
end
if A(NROW(n),n)==0;
disp('no solution');
return;
end
x=zeros(n,1);
x(n)=A(NROW(n),n+1)./A(NROW(n),n);
for k=n-1:-1:1;
SUM=0;
for l=k+1:n;
SUM=SUM+A(NROW(k),l)*x(l);
end
x(k)=(A(NROW(k),n+1)-SUM)./A(NROW(k),k);
end
x
Exercise

1. Construct a program which reduces a given matrix into echelon form
(i.e. upper triangular matrix) by using Gaussian elimination.
•   n=3;
•   A=[1 2 3 2;2 1 0 3;3 1 1 4];
•   for k=1:n-1;
•      for l=k+1:n;
•         if A(k,k)==0;
•            disp('there is zero');
•            break;
•         end
•         m=A(l,k)/A(k,k);
•         T=A(l,:)-m*A(k,:);
•         A(l,:)=T;
•      end
•   end
•   A
Exercise

2. For an nxn matrix A, construct a program such that for any specified
kth column, find the entry in this column with the biggest absolute value,
and switch the row containing this entry with the first row of A.
Exercise

3. For an nx(n+1) matrix A of echelon form which
corresponds to a system of linear equations with n
unknowns, construct a program to find the
unknowns.
Jacobi Iteration

To solve Ax=b given an initial approximation x0:
Input the number of equations and unknowns n; the entries aij,
1<=i,j<=n of the matrix A; the entries bi, 1<=i<=n; the entries
XOi , 1<=i<=n of X); tolerance TOL; maximum number of
iterations N.
Output the approximate solutions x1,x2,…,xn or a message that
the number of iterations was exceeded.
Step 1 Set k=1
Step 2 While (k<=N) do Steps 3—6
Step 3 For i=1,…,n
set

Step 4 If ||x-XO||<TOL the OUTPUT(x1,x2,…,xn);
(successful)
Stop;
Step 5 Set k=k+1;
Step 6 For i=1,2,…,n set XOi=xi.
Step 7 OUTPUT(`Maximum number of iterations exceeded’)
Gauss-Seidel Iterative

To solve Ax=b given an initial approximation x0:
Input the number of equations and unknowns n; the entries aij,
1<=i,j<=n of the matrix A; the entries bi, 1<=i<=n; the entries
XOi , 1<=i<=n of X); tolerance TOL; maximum number of
iterations N.
Output the approximate solutions x1,x2,…,xn or a message that
the number of iterations was exceeded.
Step 1 Set k=1
Step 2 While (k<=N) do Steps 3—6
Step 3 For i=1,…,n
set

Step 4 If ||x-XO||<TOL the OUTPUT(x1,x2,…,xn);
(successful)
Stop;
Step 5 Set k=k+1;
Step 6 For i=1,2,…,n set XOi=xi.
Step 7 OUTPUT(`Maximum number of iterations exceeded’)
STOP.
A=[1 0 0;0 1 0;0 0 1];
b=[2,3,4];
n=3;
XO=[1,2,3];
TOL=1e-4;
N=50;
k=1;
while k<=N;
for i=1:n;
SUM1=0;
for j=1:i-1;
SUM1=SUM1+A(i,j)*x(j);
end
SUM2=0;
for j=i+1:n;
SUM2=SUM2+A(i,j)*XO(j);
end
x(i)=1/A(i,i)*(-SUM1-SUM2+b(i));
end
if norm(x-XO)<TOL;
x
break;
end
k=k+1;
for i=1:n;
XO(i)=x(i);
end
end
if k>N;
disp('Maximum number of iterations exceeded');
end
Excercise

Construct a function GaussSeidel such that
GaussSeidel(A,b)
yields the solution to the system of linear equations Ax=b by making use of the
Gauss-Seidel iterative.
SOR

To solve Ax=b given r an initial approximation x0:
Input the number of equations and unknowns n; the entries aij,
1<=i,j<=n of the matrix A; the entries bi, 1<=i<=n; the entries
XOi , 1<=i<=n of X); the parameter r; tolerance TOL;
maximum number of iterations N.
Output the approximate solutions x1,x2,…,xn or a message that
the number of iterations was exceeded.
Step 1 Set k=1
Step 2 While (k<=N) do Steps 3—6
Step 3 For i=1,…,n
set

Step 4 If ||x-XO||<TOL the OUTPUT(x1,x2,…,xn);
(successful)
Stop;
Step 5 Set k=k+1;
Step 6 For i=1,2,…,n set XOi=xi.
Step 7 OUTPUT(`Maximum number of iterations exceeded’)
STOP.

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 3 posted: 10/19/2011 language: English pages: 32