Mетод исключения Гаусса без выбора гл. элемента на Matlab

[Alenka_kiss]

Мне нужно на языке системы MATLAB реализовать метод исключения Гаусса без выбора главного элемента в виде процедуры.
Я решала СЛАУ методом Гаусса, но обычным способом (в тетрадке) и, правда, с выбором главного элемента.
В MATLABе я нашала способ сделать расчеты с помощью стандартного решателя системы x=A\b
Вот пример. Но это, конечно, я понимаю, не то, что нужно.

» A=[1 2 3; 5 4 3; 3 4 3];
» b=[3; 9; 7];
» x=A\b

x =

    1.0000
    1.0000
    0.0000

А как составить процедуру на языке MATLAB? Я примерно представляю, как это сделать на Delphi, есть алгоритм, а вот в MATLABе...

[photon]

Alenka_kiss, брр. В чем проблема? - в алгоритме? в том как записать в MatLAB? или в том, как в MatLAB оформляются процедуры/функции?

Добавлено спустя 1 минуту 12 секунд:

Если есть реализация на Delphi, выложите тут, а мы попробуем помочь с переводом под MatLAB

[Alenka_kiss]

Проблема в том как записать в Matlabe. Алгоритм есть, выкладываю. Реализации на Delphi нет, просто я думаю ,что примерно знаю, как этот алгоритм реализовать на Delphi.
Алгоритм_рисунок

[photon]

извините, у меня сейчас нет времени разбираться в методе - я не знаю, что, например, подразумевает "перестановка уравнений". Если Вы напишете конкретную реализацию на Delphi, C, Java или чем-то что хотя бы немного на это похожем, то я могу постараться переложить на MatLAB

[Alenka_kiss]

Вот есть на С++.

Gauss.cpp - наверное, что-то похожее на это будет.

#include<iostream.h>
#include<math.h>

const n=4;
int IOR[n],k,l,i,p,M,j; 
double AKK,AMAIN,x[n],delitel;

void main()
{
double array[n][n]={4.11,-1.26,-5.99,1.29,-1.26,2.00,4.00,0.00,3.18,-1.97,0.49,-1.00,1.29,3.81,-1.56,0.00};
double vek[n]={-0.75,1.08,3.38,0.87};
/*присвоение компонентам массива перестановок IOR(k) исходных значений*/
for (k=0; k<n; k++)
{
IOR[k]=k;
}
/*нахождение индекса p*/
for (k=0; k<n; k++)
{
AKK=0;
for (i=k; i<n; i++)
{
l=IOR[i];
if (fabs(array[l][k])<AKK) continue;
M=l; p=i; AKK=fabs(array[l][k]);
}
/*меняем местами значения IOR[p]  и IOR[k], если IOR[p] не равно IOR[k]*/
AMAIN=array[M][k];
if (k!=p) {IOR[p]=IOR[k]; IOR[k]=M;};
if (AMAIN==0) {cout<<endl<<"Oshibka IER=1, programma bila avariino zavershena"; break;};
/*Исключение переменной Xk (прямой ход метода гаусса)*/
for (j=k; j<n; j++)
{
array[M][j]=array[M][j]/AMAIN;
}
vek[M]=vek[M]/AMAIN;

for (i=k+1; i<n; i++)
{
l=IOR[i];
delitel=array[l][k];
for (j=k; j<n; j++)
{
array[l][j]=array[l][j]-delitel*array[M][j];
}
vek[l]=vek[l]-delitel*vek[M]; 
}
if (array[l][n-1]==0) {cout<<endl<<"Oshibka IER=2, programma bila avariino zavershena"; break;};
}
/*Обратный ход метода Гаусса*/
double sum;
for (k=n-1; k>=0; k--)
{
l=IOR[k];
sum=0;
for (j=k+1; j<n; j++)
{
sum=sum+array[l][j]*x[j];
}
x[k]=vek[l]-sum;
}
cout<<endl<<"Vektor reshenia: ";
for (i=0; i<n; i++)
{
cout<<endl<<x[i]<<"     "<<endl;
}
/*Вычисление вектора невязки*/
vek[0]=-0.75-(x[0]*4.11-x[1]*1.26-x[2]*5.99+x[3]*1.29);
vek[1]=1.08-(-x[0]*1.26+x[1]*2.00+x[2]*4.00);
vek[2]=3.38-(x[0]*3.18-x[1]*1.97+x[2]*0.49-x[3]);
vek[3]=0.87-(x[0]*1.29+x[1]*3.81-x[2]*1.56);
cout<<endl<<endl<<"Vektor neviazki : "<<endl;
AKK=0;
for (i=0; i<n; i++)
{
if (fabs(vek[i]>AKK)) AKK=fabs(vek[i]);
cout<<vek[i]<<endl;
}
cout<<endl<<endl<<"Norma vektora neviazki : "<<AKK<<endl;
}

[photon]

два вопроса:
1) что делает fabs() ?
2) что Вы получаете на выходе для указанных данных - чтобы я мог сравнить результат.

[Alenka_kiss]

после выполнения программы получается:
Vektor reshenia:

0.799519

0.142201

0.450748

-0.896799

Vektor neviazki:
3.957338е-16
8.261621e-17
9.207045е-16
3.350727е-16

Norma vektora neviazki: 9.207045е-16

На Matlab мне нужно только решение получить (ответ).
Данные ввести при помощи генератора случайных чисел. Процедуру нужно оттестировать на СЛАУ с хорошо обусловленной матрицей 50х50 и СЛАУ с плохообусловленной матрицей 50х50. Результат сравнить с результатами расчетов с помощью стандартного решателя системы. Думаю, я теперь понятней объяснила.

На Matlab.

clc;

n = input('Введите размерность матрицы n = ');
clc;

disp('Выбирите параметры генерации матрицы:')
disp('  1 - хорошо обусловленная матрица коэффициентов')
disp('  2 - плохо обусловленная матрица коэффициентов')

if input('-> ') == 1,
  A = rand(n);
else  
  A = hilb(n);
end;

B = randn(n, 1);
A1 = A;
B1 = B;
k = 1;

while k < n, % начало 1 цикла
  ii = k + 1;
  if (A(k, k) == 0) & (k < n), % переставляем уравнения если элемент A(k, k) равен 0
    d = A(k, :);
    A(k, :) = A(k + 1, :);
    A(k + 1, :) = d;
    d = B(k);
    B(k) = B(k + 1);
    B(k + 1) = d;
  end; % конец перестановки уравнений
  
  while ii <= n, % второй цикл
    c = A(ii, k) / A(k, k);
    A(ii, k) = 0;
    jj = k + 1;
    while jj <= n, % начало третьего цикла
      A(ii, jj) = A(ii, jj) - c * A(k, jj);
      jj = jj + 1;
    end; % Конец третьего цикла

    B(ii) = B(ii) - c * B(k);
    ii = ii + 1;
  end; % конец второго цикла

  k = k + 1;
end; % конец 1 цикла

% Конец прямого хода, начинаем обратный ход

X(n) = B(n) / A(n, n);
ii = n - 1;

while ii >= 1,
  jj = ii + 1;
  S = 0;
  
  while jj <= n,
    S = S + A(ii, jj) * X(jj);
    jj = jj + 1;
  end;

  X(ii) = (B(ii) - S) / A(ii, ii);
  ii = ii - 1;
end;
clc;
A = A1;
B = B1;
A
B
disp('Результат решение уравнений с помощью A\B')
 A\B
disp('Результат решения уравнений с помощью метода Гаусса')
X'

[_Dmitry_]

# Прямой ход
for k=1:n
aa=B(k,k);
B(k,k:n+1)=B(k,k:n+1)/aa;
for i=k+1:n
   bb=B(i,k);
   B(i,k:n+1)=B(i,k:n+1)-bb*B(k,k:n+1);
end
end

# Обратный ход
x(n,1)=B(n,n+1);
for k=n-1:-1:1
  x(k,1)=(B(k,n+1)-B(k,k+1:n)*x(k+1:n,1));
end