|
Текст подпрограммы и версий zf31r_c.zip zf31d_c.zip |
Тексты тестовых примеров tzf31r_c.zip tzf31d_c.zip |
Решение системы n нелинейных уравнений с n неизвестными методом Ньютона.
Пусть дана система N нелинейных уравнений с N неизвестными.
fi ( x1, x2, ..., xN ) = 0 , i = 1, 2, ..., N
Для решения заданной системы подпрограмма zf31r_c на k - м шаге итерационного процесса Ньютона вычисляет решение следующей линейной системы методом Гаусса с выбором ведущего элемента по столбцу:
N
∑ αi j δx j( k ) = β i ,
j =1
где αi j = ∂fi /∂xj - якобиан системы, и βi = - fi, а затем вычисляет (k + 1) - е приближение к вектору решения на основе ранее найденного k - го приближения:
x ik+1 = x i(k) + δx i(k) , i = 1, 2, ..., N
Итерационный процесс по k повторяется до тех пор, пока не будет выполнен один из двух критериев сходимости:
N
(1) ∑ | δx i(k) | ≤ EPSX
i =1
N
(2) ∑ | f i ( x1(k+1), x2(k+1), ..., xN(k+1)) | ≤ EPSF
i =1
Если же ни одно из этих двух условий не выполнено за ITMAX итераций, то zf31r_c прекращает свою работу и выдает диагностическое сообщение
Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1979.
int zf31r_c (S_fp fun, real *root, integer *n, real *epsx,
real *epsf, integer *itmax, real *alpha, real *beta, integer *nlead,
integer *ierr)
Параметры
| fun - |
имя подпрограммы вычисления
fi (x1, x2,
..., xn) и ∂fi /∂xj ǀ (x1, x2, ..., xn) , i=1,2,...,n, j=1,2,...,n; оформляется в виде подпрограммы с тремя формальными параметрами: |
| x - | вещественный вектор длины n, содержащий текущее приближение к решению исходной системы; |
| alpha - |
двумерный вещественный вектор размеров n на
n , в который помещаются значения элементов
Якобиана системы в точке x, т.е. ∂fi /∂xj ǀ (x1, x2, ..., xn) ; |
| beta - | вещественный вектор длины n, в который помещаются значения - f i ( x1, x2, ..., xn ) |
| root - | вещественный вектор длины n, содержащий вычисленное решение системы; перед началом работы подпрограммы вектор root должен содержать начальное приближение к решению; |
| n - | заданное число уравнений системы (тип: целый); |
| epsx - | первый критерий сходимости (1) (тип: вещественный); |
| epsf - | второй критерий сходимости (2) (тип: вещественный); |
| itmax - | целая переменная, значение которой перед началом работы подпрограммы должно быть положено равным максимальному числу итераций, ориентировочно требуемых для обеспечения сходимости в соответствии с заданными критериями; |
| alpha - | двумерный вещественный массив размеров n на n, используемый в качестве рабочего; |
| beta - | вещественный вектор длины n, используемый в качестве рабочего; |
| nlead - | целый вектор длины n, используемый в качестве рабочего; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженный в ходе работы подпрограммы; при этом: |
| ierr=65 - | когда решение системы не может быть найдено в пределах заданного числа итераций; |
| ierr=66 - | когда якобиан заданной системы имеет особенность. |
Версии
| zf31d_c - | решение системы n нелинейных уравнений с n неизвестными методом Ньютона в режиме удвоенной точности; при этом параметры root, epsx, epsf, alpha, beta и формальные параметры подпрограммы fun должны иметь тип double. |
Вызываемые подпрограммы
|
asg9r_c - asg9d_c | решение системы Ax = b или ATx = b методом Гаусса с выбором ведущего элемента по столбцу; используются в подпрограммах zf31r_c и zf31d_c соответственно; |
|
utzf10_c - utzf11_c | подпрограммы выдачи диагностических сообщений при работе подпрограмм zf31r_c и zf31d_c соответственно. |
Замечания по использованию
| В подпрограммах zf31r_c и zf31d_c имеется внешняя структура с именем zf31rr_ , содержащая элемент целого типа с именем iter. Переменная iter полагается равной количеству итераций, выполненных при решении системы. |
struct {
int iter;
} zf31rr_;
#define zf31rr_1 zf31rr_
int main(void)
{
/* Initialized data */
static float root[3] = { 1.f,1.f,1.f };
/* Local variables */
static float beta[3], epsf;
static int ierr;
extern int zf31r_c(U_fp, float *, int *, float *, float *, int *,
float *, float *, int *, int *);
static float epsx;
static int n, nlead[3];
static float alpha[9] /* was [3][3] */;
static int itmax;
extern int fun_c();
n = 3;
epsx = 1e-5f;
epsf = 1e-5f;
itmax = 30;
zf31r_c((U_fp)fun_c, root, &n, &epsx, &epsf, &itmax, alpha, beta, nlead,
&ierr);
printf("\n %15.6e %15.6e %15.6e \n", root[0], root[1], root[2]);
printf("\n %5i \n", zf31rr_1.iter);
printf("\n %5i \n", ierr);
return 0;
} /* main */
int fun_c(float *x, float *alpha, float *beta)
{
/* System generated locals */
float r__1;
/* Builtin functions */
double exp(double), sin(double), cos(double);
#define alpha_ref(a_1,a_2) alpha[(a_2)*3 + a_1]
/* Parameter adjustments */
--beta;
alpha -= 4;
--x;
/* Function Body */
/* Computing 2nd power */
r__1 = x[2] + x[3];
beta[1] = -(x[1] + (float)exp(x[1] - 1) + r__1 * r__1 - 27);
/* Computing 2nd power */
r__1 = x[3];
beta[2] = -(x[1] * (float)exp(x[2] - 2) + r__1 * r__1 - 10);
/* Computing 2nd power */
r__1 = x[2];
beta[3] = -(x[3] + (float)sin(x[2] - 2) + r__1 * r__1 - 7);
alpha_ref(1, 1) = (float)exp(x[1] - 1) + 1;
alpha_ref(1, 2) = (x[2] + x[3]) * 2;
alpha_ref(1, 3) = (x[2] + x[3]) * 2;
alpha_ref(2, 1) = (float)exp(x[2] - 2);
alpha_ref(2, 2) = x[1] * (float)exp(x[2] - 2);
alpha_ref(2, 3) = x[3] * 2;
alpha_ref(3, 1) = 0.f;
alpha_ref(3, 2) = (float)cos(x[2] - 2) + x[2] * 2;
alpha_ref(3, 3) = 1.f;
return 0;
} /* fun_c */
Результаты: root = ( 1.0, 2.0, 3.0 )
iter = 7 , ierr = 0