echo on; clc;
%---------------------------------------------------------------------------
%A2_9 MATLAB script file for implementing Algorithm 2.9
%
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1995
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
% Prentice Hall, International Editions: ISBN 0-13-625047-5
% This free software is compliments of the author.
% E-mail address: in%"mathews@fullerton.edu"
%
% Algorithm 2.9* (Fixed Point Iteration in higher dimensions).
% Section 2.6, Iteration for Nonlinear Systems, Page 108
%---------------------------------------------------------------------------
clc; clear all; format long;
% - - - - - - - - - - - - - - - - - - - - - - -
%
% This program implements fixed point iteration
% in 2-dimensions for solving X = G(X).
% Use the vector notation X = (x y).
%
% Define and store the functions g1(X) and g2(X) and
% define and store the function G(X) in the
% M-files g1.m g2.m and G.m respectively.
%
% function z = g1(x,y)
% z = (x.^2 - y + 0.5)./2;
% function z = g2(x,y)
% z = (- x.^2 - 4.*y.^2 + 8.*y + 4)./8;
% function Z = G(X)
% x = X(1); y = X(2);
% Z = [g1(x,y) g2(x,y)];
pause % Press any key to continue.
clc;
%.......................................................................
% Begin a section which enters the function(s) necessary for the example
% into M-file(s) by executing the diary command in this script file.
% The preferred programming method is not to use these steps.
% One should enter the function(s) into the M-file(s) with an editor.
delete output
delete g1.m
diary g1.m; disp('function z = g1(x,y)');...
disp('z = (x.^2 - y + 0.5)./2;');...
diary off;
delete g2.m
diary g2.m; disp('function z = g2(x,y)');...
disp('z = (- x.^2 - 4.*y.^2 + 8.*y + 4)./8;');...
diary off;
delete G.m
diary G.m; disp('function Z = G(X)');...
disp('x = X(1); y = X(2);');...
disp('Z = [g1(x,y) g2(x,y)];');...
diary off;
% Remark. g1.m, g2.m, G.m and fix2dim.m are used for Algorithm 2.9
g1(0,0); g2(0,0); G([0 0]); % Test for files g1.m g2.m G.m
pause % Press any key to graph x = g1(x,y) and y = g2(x,y).
clc;
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
% Prepare graphics arrays to plot x = g1(x,y) and y = g2(x,y).
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
a = -2.0;
b = 2.5;
c = -1.0;
d = 1.5;
hx = (b-a)/31;
hy = (d-c)/31;
[X Y] = meshgrid(a:hx:b, c:hy:d);
Z = g1(X,Y) - X;
W = g2(X,Y) - Y;
clc; figure(1); clf;
%~~~~~~~~~~~~~~~~~~~~~~~
% Begin graphics section
%~~~~~~~~~~~~~~~~~~~~~~~
a = -2.0;
b = 2.5;
c = -1.0;
d = 1.5;
whitebg('w');
plot([a b],[0 0],'b',[0 0],[c d],'b');
axis([a b c d]);
axis(axis);
hold on;
contour(X,Y,Z,[0 0],'-g');
contour(X,Y,W,[0 0],'-g');
grid;
title('Implicit plot of x = g1(x,y) and y = g2(x,y).');
hold off;
figure(gcf); pause % Press any key to perform fixed point iteration.
clc;
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
% Example, page 101 Use fixed point iteration in 2 dimensions
% to solve X = G(X). The system of equations is
%
% x = (x^2 - y + 0.5)/2
% y = (- x^2 - 4y^2 + 8y + 4)/8
%
% Enter the starting value in [p0 q0]
% Enter the tolerance in delta
% Enter the number of iterations in max1
p0 = 0.0;
q0 = 1.0;
P0 = [p0 q0];
delta = 1e-12;
max1 = 50;
[P0,err,P] = fix2dim('G',P0,delta,max1);
pause % Press any key for the list of iterations.
clc;
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
% Prepare arrays to graph and print results.
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
a = -0.275;
b = 0.025;
c = 0.99;
d = 1.002;
hx = (b-a)/20;
hy = (d-c)/20;
[X Y] = meshgrid(a:hx:b, c:hy:d);
Z = g1(X,Y) - X;
W = g2(X,Y) - Y;
X0 = P(:,1);
Y0 = P(:,2);
clc; figure(2); clf;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Begin graphics section for the results.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
a = -0.275;
b = 0.025;
c = 0.99;
d = 1.002;
whitebg('w');
plot([a b],[0 0],'b',[0 0],[c d],'b');
axis([a b c d]);
axis(axis);
hold on;
contour(X,Y,Z,[0 0],'-g');
contour(X,Y,W,[0 0],'-g');
plot(X0,Y0,'or');
title('Graphical presentation of the fixed point iteration.');
grid;
hold off;
figure(gcf); pause % Press any key to continue.
clc;
%............................................
% Begin section to print the results.
% Diary commands are included which write all
% the results to the Matlab textfile output
%............................................
Mx1 = 'Computations for the fixed point iteration method.';
Mx2 = ' p(k) q(k)';
Mx3 = 'The solution is:';
Mx4 = 'The error estimate for P is ';
Mx5 = num2str(err);
Mx6 = [Mx4,'[~ ',Mx5,', ~ ',Mx5,']'];
clc,echo off,diary output,...
disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
disp('Iteration is converging linearly to the root.'),...
disp(''),disp(Mx3),disp(''),disp('G(P) = P = '),disp(P0),...
disp(''),disp([Mx6]),diary off,echo on
pause % Press any key to perform fixed point iteration.
clc;
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
% Example, page 101 Use fixed point iteration in 2 dimensions
% to solve X = G(X). The system of equations is
%
% x = (x^2 - y + 0.5)/2
% y = (- x^2 - 4y^2 + 8y + 4)/8
%
% Enter the starting value in [p0 q0]
% Enter the tolerance in delta
% Enter the number of iterations in max1
p0 = 2.0;
q0 = 0.0;
P0 = [p0 q0];
delta = 1e-12;
max1 = 6;
[P0,err,P] = fix2dim('G',P0,delta,max1);
pause % Press any key for the list of iterations.
clc;
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
% Prepare arrays to graph and print results.
% ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
a = -2;
b = 10;
c = -3;
d = 1;
hx = (b-a)/31;
hy = (d-c)/31;
[X Y] = meshgrid(a:hx:b, c:hy:d);
Z = g1(X,Y) - X;
W = g2(X,Y) - Y;
clc; figure(3); clf;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Begin graphics section for the results.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
a = -2;
b = 10;
c = -3;
d = 1;
whitebg('w');
plot([a b],[0 0],'b',[0 0],[c d],'b');
axis([a b c d]);
axis(axis);
hold on;
contour(X,Y,Z,[0 0],'-g');
contour(X,Y,W,[0 0],'-g');
X0 = P(:,1);
Y0 = P(:,2);
plot(X0,Y0,'or');
plot([a b],[0 0],'b',[0 0],[c d],'b');
title('Graphical presentation of the fixed point iteration.');
grid;
hold off;
figure(gcf); pause % Press any key to continue.
clc;
%............................................
% Begin section to print the results.
% Diary commands are included which write all
% the results to the Matlab textfile output
%............................................
clc,echo off, diary output,...
disp(''), disp(Mx1),disp(''), disp(Mx2),disp(P),...
disp('Iteration is diverging to "infinity".'),...
disp('Note the "scale factor" for the table of computations.'),...
diary off, echo on
