9 views (last 30 days)

The figure below shows Newtons Raphson Iteration to achieve the desired accuracy along x axis using the intensity of blue color tone . How can I make such figure I have the following code for achieving the basin of attraction. Please also guide that how can I find that how many iterations the program is doing to reach a certain attractor for each initial condition.

clc ; clear all

warning('off') % To off the warning which shows "Matrix is close to singular

% badly scaled" when algorithm passes through a point where the Jacobian

% matrix is singular

% The roots of the given governing equations

r1 = [-1 ;-1] ;

r2 = [0 ;0] ;

r3 = [1 ;1] ;

% Initial conditions

x = linspace(-2,2,200) ;

y = linspace(-2,2,200) ;

% Initialize the required matrices

Xr1 = [] ; Xr2 = [] ; Xr3 = [] ; Xr4 = [] ;

tic

for i = 1:length(x)

for j = 1:length(y)

X0 = [x(i);y(j)] ;

% Solve the system of Equations using Newton's Method

X = NewtonRaphson(X0) ;

% Locating the initial conditions according to error

if norm(X-r1)<1e-8

Xr1 = [X0 Xr1] ;

elseif norm(X-r2)<1e-8

Xr2 = [X0 Xr2] ;

elseif norm(X-r3)<1e-8

Xr3 = [X0 Xr3] ;

else % if not close to any of the roots

Xr4 = [X0 Xr4] ;

end

end

end

toc

warning('on') % Remove the warning off constraint

% Initialize figure

figure

set(gcf,'color','w')

hold on

plot(Xr1(1,:),Xr1(2,:),'.','color','r') ;

plot(Xr2(1,:),Xr2(2,:),'.','color','b') ;

plot(Xr3(1,:),Xr3(2,:),'.','color','g') ;

plot(Xr4(1,:),Xr4(2,:),'.','color','k') ;

title('Basin of attraction for f(x,y) = x^3-y = 0 and y^3-x=0')

and the Newton Raphson Method file is

function X = NewtonRaphson(X)

NoIter = 10 ;

% Run a loop for given number of iterations

for j=1:NoIter

% Governing equations

f=[X(1)^3 - X(2); X(2)^3 - X(1)] ;

% Jacobian Matrix

Jf=[3*X(1)^2 -1; -1 3*X(2)^2];

% Updating the root

X=X-Jf\f;

end

J. Alex Lee
on 21 Aug 2020

Organize your "list" of initial guesses (x,y) as a meshgrid

[X,Y] = meshgrid(x,y)

Then easier to loop through the "unwrapped" arrays

for k = 1:numel(X)

X0 = [X(k);Y(k)];

end

Then you can plot any output of your solution as an "image".

Also, not sure your current code actually does what you are asking...you don't seem to be recording the number of iterations, and it doesn't seem to me that you need to know which roots you are converging to in order to create the first plot...

J. Alex Lee
on 24 Aug 2020

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.