Решение СНАУ методом Ньютона

[Tina]

если у кого-нибудь есть исходник ~решения системы нелинейных уравнений методом Ньютона~, пожалуйста подскажите
у меня он есть на С++, но хотелось бы на Паскале,т.к. в С++ я не очень понимаю



спасибо

[ZhekON]

так скинь на С++ а знатоки переведут))

[Cannibal]

вот исходник. подставь свою функцию

program hord;
uses crt;
var b,x1,x2,e:real;

function f(x0:real):real;
Var g:real;
begin g:={Функция};
      F:=g
end;
function f1(x0:real):real;
var  g1:real;
begin g1:={производная};
      f1:=g1;
end;

begin Write('Введите первое приближение '); read(x1);
      write('Введите неподвижную точку '); read(b);
      write('введите точность '); read (e);
      repeat
            x2:=x1;
            X1:=x1-F(x1)/F1(x1);
            writeln (x1);
      until abs(x1-x2)<e;
      writeln (X1);
end.

[_Dmitry]

const
  n_max = 2; //число уравнений

type
  TVec = array[1..n_max] of double;
  TMat = array[1..n_max,1..n_max] of double;

//система уравнений
procedure F(n: integer; x: TVec; var y: TVec);
begin
  y[1]:=x[1]-exp(-x[2]);
  y[2]:=x[2]-exp(x[1]);
end;

//вычисление обратной матрицы Якоби
procedure G(n: integer; x: TVec; var a: TMat);
var
  det: double;
  i,j: integer;
begin
  a[1,1]:=1; a[1,2]:=-exp(-x[2]);
  a[2,1]:=exp(x[1]); a[2,2]:=1;
  det:=1+exp(x[1])*exp(-x[2]);
  for i:=1 to n do
    for j:=1 to n do
      a[i,j]:=a[i,j]/det;
end;

//процедура решения системы методом Ньютона
procedure Newts(n: integer; var x: TVec; eps: double; var k: integer);
var
  i,j: integer;
  x0,y: TVec;
  a: TMat;
  cod: boolean;
begin
  k:=0;
  repeat
    for i:=1 to n do x0[i]:=x[i];
    F(n,x,y); G(n,x,a);
    for i:=1 to n do
      for j:=1 to n do
        x[i]:=x[i]-a[i,j]*y[j];
    cod:=true;
    for i:=1 to n do
      if abs(x[i]-x0[i]) > eps then cod:=false;
    k:=k+1;
  until cod;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  x: TVec;
  n,k: integer;
  eps: double;
begin
  n:=n_max; //число уравнений
  eps:=0.001; //точность
  x[1]:=0.5; x[2]:=0.5; //начальное приближение
  Newts(n,x,eps,k); //решение системы
  Memo1.Lines.Add('Решение системы:');  
  Memo1.Lines.Add('x[1]='+FloatToStr(x[1]));
  Memo1.Lines.Add('x[2]='+FloatToStr(x[2]));
  Memo1.Lines.Add('Число итераций: '+IntToStr(k));  
end;

Естественно, процедуры F и G нужно написать свои.

[Tina]

Всем огромное спасибо)))))

[Tina]

Всем огромное спасибо)))))

[MrVenom]

Всем доброе время суток.
Вот требуется решить такую систему:
1.5160513 + 1.6 * x1 + 2.2 * Cos(2.2 * x2) + 1.6 * Sin(1.1 * x3) - 1.6 * x4^4=0
2.9679284 + 1.8 * x2^3 - 2.2 * Cos(0.8 * x2) - 2.8 * Exp(2.5*x2)=0
0.9979576 + 2.1 * x2^3 - 2.5 * x3 + 2.3 * x3^2 + 2 * x4^2 + 0.2 * Sin(-0.4 * x4)=0
-0.7499998 + 2.5 * x1^2 + 1.7 * x1^3 + 0.6 * x3 + 1.2 * x4^3=0
Начальные условия: x1 = -1.27; x2 = -0.71; x3 = 0.13; x4 = -0.7;
методом Ньютона.
Вот что наваял, причем для уравнения(которое в комментах все работает норм, а если взять свое, то нифига )
К сожалению, код во вложении полный, т.к. больше 5000 символов, тут выложу без определения детерминанта матрицы:

        static void F(double[] x, ref double[] y)//уравнения
        {
            y[0] = 1.5160513 + 1.6 * x[0] + 2.2 * Math.Cos(2.2 * x[1]) + 1.6 * Math.Sin(1.1 * x[2]) - 1.6 * Math.Pow(x[3], 4);
            y[1] = 2.9679284 + 1.8 * Math.Pow(x[1], 3) - 2.2 * Math.Cos(0.8 * x[1]) - 2.8 * Math.Exp(2.5*x[1]);
            y[2] = 0.9979576 + 2.1 * Math.Pow(x[1], 3) - 2.5 * x[2] + 2.3 * Math.Pow(x[2], 2) + 2 * Math.Pow(x[3], 2) + 0.2 * Math.Sin(-0.4 * x[3]);
            y[3] = -0.7499998 + 2.5 * Math.Pow(x[0], 2) + 1.7 * Math.Pow(x[0], 3) + 0.6 * x[2] + 1.2 * Math.Pow(x[3], 3);
            /*
            y[0] = x[0] - Math.Exp(-x[1]);
            y[1] = x[1] - Math.Exp(x[0]);*/
            }
        static void Jacobian(double[] x, ref double[][] a)
        {
            double det;
            //забиваем матрицу частными производными            
            a[0][0] = 1.6; a[0][1] = -(121 * Math.Sin(2.2 * x[1]))/25;
            a[0][2] = (44 * Math.Cos(1.1 * x[2])) / 25; a[0][3] = -6.4 * Math.Pow(x[3],3);
            a[1][0] = 0; a[1][1] = -7 * Math.Exp(2.5 * x[1]) + (44 * Math.Sin(0.8 * x[1])) / 25 + 5.4 * Math.Pow(x[1], 2);
            a[1][2] = 0; a[1][3] = 0;
            a[2][0] = 0; a[2][1] = 6.3 * Math.Pow(x[1], 2); a[2][2] = 4.6 * x[2] - 2.5;
            a[2][3] = 4 * x[3] - (2 * Math.Cos(0.4 * x[3])) / 25;
            a[3][0] = 5.1 * Math.Pow(x[0], 2) + 5 * x[0]; a[3][1] = 0; a[3][2] = 0.6; a[3][3] = 3.6 * Math.Pow(x[3], 2);
            det = Det(a, 4);
            for (int i = 0; i < 4; i++)
                for (int j = 0; j < 4; j++) a[i][j] = a[i][j] / det;
            /*
            double det;
            a[0][0]=1; a[0][1]=Math.Exp(-x[1]);
            a[1][0]=Math.Exp(x[0]); a[1][1]=1;
            det=1+Math.Exp(x[0])*Math.Exp(-x[1]);
            for (int i = 0; i < 2; i++)
                for (int j = 0; j < 2; j++) a[i][j] = a[i][j] / det;*/
        }

        static void Newthon(ref double[] x, double eps, ref int k)
        {
            bool cod;
            double[] x0 = new double[4];
            double[] y = new double[4];
            double[][] a = new double[4][];
            for (int i = 0; i < 4; i++) a[i] = new double[4];
            do
            {
                for (int i = 0; i < 4; i++) x0[i] = x[i];
                F(x, ref y);
                Jacobian(x, ref a);
                for (int i = 0; i < 4; i++)
                    for (int j = 0; j < 4; j++) x[i] = x[i] - a[i][j] * y[j];
                cod = true;
                for (int i = 0; i < 4; i++)
                    if (Math.Abs(x[i] - x0[i]) > eps) cod = false;
                k++;
            } while (cod == false);
            
            /*
            for (int i = 0; i < 2; i++) a[i] = new double[2];
            do
            {
                for (int i = 0; i < 2; i++) x0[i] = x[i];
                F(x, ref y); Jacobian(x, ref a);
                for (int i = 0; i < 2; i++)
                    for (int j = 0; j < 2; j++) x[i] = x[i] - a[i][j] * y[j];
                cod = true;
                for (int i = 0; i < 2; i++)
                    if (Math.Abs(x[i] - x0[i]) > eps) cod = false;
                k++;
            } while (cod == false);*/
        }
        static void Main(string[] args)
        {
            double eps = 0.001;
            double[] x = new double[4];
            x[0] = -1.27; x[1] = -0.71; x[2] = 0.13; x[3] = -0.7;
            int k = 0;
            Newthon(ref x, eps, ref k);
            for (int i = 0; i < 4; i++) Console.WriteLine(x[i]);//выводим иксы...
            Console.WriteLine(k);
            double[] y = new double[4];
            F(x, ref y);
            for (int i = 0; i < 4; i++) Console.WriteLine(y[i]);//это для проверки...
            /*
            double eps = 0.00001;
            double[] x = new double[2];
            x[0] = 0.5; x[1] = 0.5;
            int k = 0;
            Newthon(ref x, eps, ref k);
            for (int i = 0; i < 2; i++) Console.WriteLine(x[i]);
            Console.WriteLine(k);
            double[] y = new double[2];
            F(x, ref y);
            for (int i = 0; i < 2; i++) Console.WriteLine(y[i]);*/
            Console.ReadKey();
        }

Вопрос решен

[mrsikilinda]

Доброе утро. У меня возникла проблема с решением СНАУ. Есть код Метода Ньютона необходимо написать код модифицированным методом Ньютона на Delphi. Вы не могли бы написать код модифицированного метода Ньютона для примера Дмитрия. Заранее благодарен