|
ru.algorithms
- RU.ALGORITHMS ----------------------------------------------------------------
From : Stanislav Shwartsman 2:400/520 07 Feb 2002 20:12:57
To : Eduard Vatutin
Subject : Решение уравнения
--------------------------------------------------------------------------------
06 Feb 02 01:29, you wrote to Evgenij Masherov:
EV> Можно хоть один из методов расписать. А то что-то мои поиски пока
EV> результата не дали. Метод Hьютона
EV> Xn+1 = Xn - f(Xn)/f'(Xn)
EV> насколько мне известно, зависит от начального приближения и сходится к
EV> одному из корней (если вообще сходится), а нужны хотя бы _все_
EV> действительные.
Метод Bairstow, находит все корни (в том числе и комплексные) с заданной
точностью. Исходник для матлаба, кто хочет - переделывайте на C/C++ ...
Писался как домашка по курсу "Численный Анализ 1".
=== Cut ===
function [u,v,b,i] = bairstaw(AVec, u0, v0, tol, NMax)
%
% Input: AVec - Vector of coefficients of polynomial P(x)
% u0,v0 - Initial approximation
% tol - Reasonable error
% Nmax - Maximum number of iterations
%
% Output: u,v - The final u and v
% b - Q(x) = P(x) / (x^2 - u*x - v)
% i - Number of iterations done
%
u = u0;
v = v0;
i = 1; % counter of iterations
n = length(AVec);
b(1) = 1; %
b(2) = 1; % to make the while happy
Px = [1 1]; %
while (i < NMax) & (norm(b) > tol) & (norm(Px) > tol)
i=i+1;
% The calculation of vector b
b = zeros(1, n+2); % by the algo b(n+1) = b(n+2) = 0
for(j=n:-1:1)
b(j) = AVec(j)+u*b(j+1)+v*b(j+2);
end;
% The calculation of vector c
c = zeros(1, n+1); % by the algo c(n) = c(n+1) = 0
for(j=n-1:-1:1)
c(j) = b(j+1)+u*c(j+1)+v*c(j+2);
end;
% J calculation
J = [c(1) c(2); c(2) c(3)];
if(abs(det(J)) < 10^-15)
disp('The J matrix is singular!');
i=NMax+1;
else
dV = -inv(J)*[b(1) b(2)]';
u = u+dV(1);
v = v+dV(2);
end;
x0(1) = (u+(u^2+4*v)^0.5)/2;
x0(2) = (u-(u^2+4*v)^0.5)/2;
Px = [AVec(1) AVec(1)];
for (k=1:n-1)
Px = Px+AVec(k+1)*[x0(1) x0(2)].^k;
end;
end;
if(i==NMax)
disp('Maximum number of iterations reached')
end
b = b(:,[3:end]);
=== Cut ===
=== Cut ===
function r=find_roots(A)
%
% This function uses Newton-Bairstaw algorithm for calculating
% roots of the given polynomial.
%
% The polynome is: p(x) = an*x^n + ... + a0
%
% The coefficients of the polynomial must be given in vector
% A = [A0, A1, A2, A3 ... An] then P(x) = An*x^n + ... + A0
%
% Output of the function is vector of the roots of P(x)
%
k=1;
while(length(A) > 2)
x = quad_root(A(:,[1 2 3])); % calc roots of sq. equation
z1=x(1);
z2=x(2);
% call teh bairstaw
[u,v,q,i] = bairstaw(A, z1+z2, -z1*z2, 10^-5, 100);
if (i>=25) % if bairstaw failed ...
x = quad_root(A(:,[3 2 1]));
z1 = 1/x(1);
z2 = 1/x(2);
[u,v,q,i] = bairstaw(A, z1+z2, -z1*z2, 10^-5, 25);
if (i>=25) % if bairstaw failed ...
disp('The method failed to find a solutions of the polynomial')
return;
end;
end;
x = quad_root([1 -u -v]);
r(k) = x(1);
r(k+1) = x(2);
k=k+2;
A = q(:,[1:end-2]); % going to the next iteration
end
if(length(A) > 1) % linear part are reminded
r(k) = -A(1)/A(2);
end;
=== Cut ===
=== Cut ===
function X=quad_root(A)
D = A(2)^2-4*A(1)*A(3);
X(1) = (-A(2)+D^0.5)/(2*A(1));
X(2) = (-A(2)-D^0.5)/(2*A(1));
=== Cut ===
E-mail: gate@fidonet.org.il
Voice Phones: 972-4-8330554 (home), 972-5-4481073 (cell)
Bye !
Stanislav (AKA Night's Man) [Team Technion]
---
* Origin: Gate From Another World ... From Haifa, Israel (2:400/520)
Вернуться к списку тем, сортированных по: возрастание даты уменьшение даты тема автор
|
