Numerical Analysis Lab

Document Sample
Numerical Analysis Lab Powered By Docstoc
					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