Basins of attraction and Newtons Method.
7 次查看(过去 30 天)
显示 更早的评论
I am working on basins of attraction for system of two equations and I need to know what is the problem of the code to get error?
the equations:
-7.716*10^-8 + 1.*10^-6*X(1)^2 +X(1)^2*(-0.01 + 0.05*X(2)^2)^2
-4.37*10^-7 + 2.14*10^-7*X(2)^2 + X(2)^2*(-0.007 + 0.21*X(1)^2 - 0.05*X(2)^2)^2
and use the code with newton method but it gives error saying'Error in BasinofAttractiontry1_3_fixedpoint (line 20)
X = NewtonRaphson(X0) ;
I don't know what is the problem exactly.
this is the code of bassin of attraction
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 = [-0.266794 ;-0.433493] ;
r2 = [0.0330718 ;0.0899008] ;
r3 = [0.260471 ;0.416437] ;
% 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') ;
the newton raphson code
function X = NewtonRaphson(X)
NoIter = 10 ;
% Run a loop for given number of iterations
for j=1:NoIter
% Governing equations
f = [-7.716*1e-8 + 1e-6*X(1)^2 + X(1)^2*(-0.01 + 0.05*X(2)^2)^2; -4.37*1e-7 + 2.14*1e-7*X(2)^2 + X(2)^2*(-0.007 + 0.21*X(1)^2 - 0.05*X(2)^2)^2];
% Jacobian Matrix
Jf=[2.*1e^-6*X(1) + 2*X(1)*(-0.01 + 0.05*X(2)^2)^2 0.197*X(1)^2*X(2)*(-0.01+ 0.05*X(2)^2); 0.846*X(1)*X(2)^2*(-0.007 + 0.21*X(1)^2 - 0.0506*X(2)^2) 4.287*1e^-7*X(2) - 0.202*X(2)^3*(-0.007 + 0.212*X(1)^2 - 0.0506*X(2)^2) + 2*X(2)*(-0.007 + 0.212*X(1)^2 - 0.0506*X(2)^2)^2];
% Updating the root
X=X-Jf\f;
end
回答(1 个)
Mischa Kim
2021-1-20
One problem was in the governing equations for f. I cleaned up a bit and this works now:
% Governing equations
f = [-7.716*1e-8 + 1e-6*X(1)^2 + X(1)^2*(-0.01 + 0.05*X(2)^2)^2;...
-4.37*1e-7 + 2.14*1e-7*X(2)^2 + X(2)^2*(-0.007 + 0.21*X(1)^2 - 0.05*X(2)^2)^2];
However, there is another problem in the plots which are probably easy to fix.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!