Уведомления

Группа в Telegram: @pythonsu

#1 Окт. 5, 2015 20:04:34

Varmint
Зарегистрирован: 2015-07-23
Сообщения: 10
Репутация: +  -1  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Пытаюсь написать скрипт, который бы решал методом Гаусса матрицы хотя бы до 10х10 включительно. Совсем не получается. Искал в Интернете готовые скрипты - не подходят даже после каких-либо изменений.

Мои попытки на 2.7.5:

from __future__ import division
matrix = [[1,1,1],
         [1,-1,2],
         [2,-1,-1]]
l = [6,
     5,
    -3]
def gauss(a,b):
 for k in range(len(a)-1):
   for i in range(k+1,len(a)):
     t = a[i][k]/a[k][k]
     for j in range(k+1,len(a)):
       a[i][j] -= a[i][k]*t
       b[i] -= b[k]*t
 # printing out results
   for i in range(len(a)):
     print(a[i],b[i])
gauss(matrix,l)

P.S. Пусть вас не смущает функция print. Я пишу на Pythonista. Там разработчик предусмотрел возможность писать print со скобочками и без.

Отредактировано Varmint (Окт. 6, 2015 20:15:29)

Офлайн

#2 Окт. 5, 2015 22:49:41

doza_and
От:
Зарегистрирован: 2010-08-15
Сообщения: 4138
Репутация: +  252  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Varmint
Искал в Интернете готовые скрипты - не подходят даже после каких-либо изменений.
Мда. А что искали?
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html



Офлайн

#3 Окт. 6, 2015 02:45:20

py.user.next
От:
Зарегистрирован: 2010-04-29
Сообщения: 9890
Репутация: +  854  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Varmint
P.S. Пусть вас не смущает функция print. Я пишу на Pythonista. Там разработчик предусмотрел возможность писать print со скобочками и без.
Вообще, через __future__ во втором питоне можно импортировать функцию print().
>>> from __future__ import print_function
>>> 
>>> print(1, 2, 3, sep='x')
1x2x3
>>>

Varmint
Пытаюсь написать скрипт, который бы решал методом Гаусса матрицы хотя бы до 10х10 включительно.
Давным-давно я делал для NxM, но так и не сделал. Сделал только приведение к треугольному виду, но система может иметь фундаментальную систему решений (когда для точных значений переменных не хватает информации).

Неплохое описание
Фундаментальная система решений

А смысл делать алгоритм для квадратной, если основная матрица системы может быть неквадратной?
Пример с вики для квадратной



Отредактировано py.user.next (Окт. 6, 2015 02:48:08)

Офлайн

#4 Окт. 6, 2015 19:42:19

Varmint
Зарегистрирован: 2015-07-23
Сообщения: 10
Репутация: +  -1  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

doza_and
Мда. А что искали?
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html
Суть как раз-таки в том, что надо самому написать, а с NumPy я не сталкивался ещё, поэтому даже сурс плохо понимаю. Если Вы можете переписать использованные в примере функции на обычный язык, то огромное Вам спасибо.
py.user.next
А смысл делать алгоритм для квадратной, если основная матрица системы может быть неквадратной?
В моём случае это только квадратная.
А ссылки эти читал, но они мало чем помогли.

UPD: не обязательно, чтобы решала 10х10.

Офлайн

#5 Окт. 6, 2015 19:45:28

Varmint
Зарегистрирован: 2015-07-23
Сообщения: 10
Репутация: +  -1  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

# Copyright (c) Isaac Evans 2011
# All rights reserved.
def myGauss(m):
    # eliminate columns
    for col in range(len(m[0])):
        for row in range(col+1, len(m)):
            r = [(rowValue * (-(m[row][col] / m[col][col]))) for rowValue in m[col]]
            m[row] = [sum(pair) for pair in zip(m[row], r)]
    # now backsolve by substitution
    ans = []
    m.reverse() # makes it easier to backsolve
    for sol in range(len(m)):
            if sol == 0:
                ans.append(m[sol][-1] / m[sol][-2])
            else:
                inner = 0
                # substitute in all known coefficients
                for x in range(sol):
                    inner += (ans[x]*m[sol][-2-x])
                # the equation is now reduced to ax + b = c form
                # solve with (c - b) / a
                ans.append((m[sol][-1]-inner)/m[sol][-sol-2])
    ans.reverse()
    return ans

Нашёл этот пример, но иногда выплывает ZeroDivisionError. А вообще, скрипт выглядит красиво.

Офлайн

#6 Окт. 6, 2015 21:40:34

doza_and
От:
Зарегистрирован: 2010-08-15
Сообщения: 4138
Репутация: +  252  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

ну на питоне факторизацию нет смысла писать.
Простейший код выглядит примерно так:

template <class TReal>
void matr_com_lux<TReal>::dcmp(void)
{
int i,j,k;
 for(k=0;k<Dim-1;k++)
 {
  for(i=k+1;i<Dim;i++)
  {
   d_alg(lu[k+k*Dim]==0.,"singularyty in dcmp");
   lu[i+k*Dim]/=lu[k+k*Dim];
   for(j=k+1;j<Dim;j++)
   {
    lu[i+j*Dim]-=lu[i+k*Dim]*lu[k+j*Dim];
   }
  }
 }
}
Но обычно обязательно требуется выделение ведущего элемента. Тогда будет чуть сложнее (это легкая модификация кода из numerical reciples)

#define TINY 1.0e-20
template <class TReal>
void matr_com_luo<TReal>::dcmp(void)
{
int i,imax,j,k;
TReal aamax,dum,sum,d,*vv;
vv=new TReal[Dim];
 d=1.;
 for(i=0;i<Dim;i++)
 {
  aamax=0.;
  for(j=0;j<Dim;j++){if(fabs(lu[i*Dim+j])>aamax)aamax=fabs(lu[i*Dim+j]);}
  d_alg(aamax==0.,"singularyty in dcmp");
  vv[i]=TReal(1.)/aamax;
 }
 for(j=0;j<Dim;j++)
 {
  for(i=0;i<j;i++)
  {
   sum=lu[i*Dim+j];
   for(k=0;k<i;k++){sum-=lu[i*Dim+k]*lu[k*Dim+j];}
   lu[i*Dim+j]=sum;
  }
  aamax=0.;
  for(i=j;i<Dim;i++)
  {
   sum=lu[i*Dim+j];
   for(k=0;k<j;k++){sum-=lu[i*Dim+k]*lu[k*Dim+j];}
   lu[i*Dim+j]=sum;
   dum=vv[i]*fabs(sum);
   if(dum>=aamax){imax=i;aamax=dum;}
  }
  if(j!=imax)
  {
   for(k=0;k<Dim;k++)
   {
    dum=lu[imax*Dim+k];
    lu[imax*Dim+k]=lu[j*Dim+k];
    lu[j*Dim+k]=dum;
   }
   d=-d;
   vv[imax]=vv[j];
  }
  indx[j]=imax;
  if(lu[j*Dim+j]==TReal(0.))lu[j*Dim+j]=TReal(TINY);
  if(j!=Dim)
  {
   dum=TReal(1.)/lu[j*Dim+j];
   for(i=j+1;i<Dim;i++){lu[i*Dim+j]=lu[i*Dim+j]*dum;}
  }
 }
delete [] vv;
}
#undef TINY
Если взять Lapack или Atlas то будет еще более кучеряво.
Если вы хотите работать с недоопределенными и переопределенными системами и получать решения то будет еще более кучеряво.



Отредактировано doza_and (Окт. 6, 2015 21:42:46)

Офлайн

#7 Окт. 7, 2015 03:15:20

py.user.next
От:
Зарегистрирован: 2010-04-29
Сообщения: 9890
Репутация: +  854  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Varmint
В моём случае это только квадратная.
А ссылки эти читал, но они мало чем помогли.

UPD: не обязательно, чтобы решала 10х10.
В том-то и дело, что если писать, то она должна решать любую систему (любые размерности с любыми результатами). А если она любую не решает, то это и не решение, а полурешение.

doza_and
for(j=0;j<Dim;j++){if(fabs(lu[i*Dim+j])>aamax)aamax=fabs(lu[i*Dim+j]);}
Где ты откопал такое?



Отредактировано py.user.next (Окт. 7, 2015 03:18:29)

Офлайн

#8 Окт. 7, 2015 04:58:11

JOHN_16
От: Россия, Петропавловск-Камчатск
Зарегистрирован: 2010-03-22
Сообщения: 3292
Репутация: +  221  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Varmint
P.S. Пусть вас не смущает функция print. Я пишу на Pythonista. Там разработчик предусмотрел возможность писать print со скобочками и без.
выкинуть это из головы, это возможности Питона начиная с какой то версии:
Python 2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)] on wi32
Type "help", "copyright", "credits" or "license" for more information.
>>> print 'qwe'
qwe
>>> print('qwe')
qwe
>>> type(print)
  File "<stdin>", line 1
    type(print)
             ^
SyntaxError: invalid syntax
>>>
>>> from __future__ import print_function
>>> print('qwe')
qwe
>>> type(print)
<type 'builtin_function_or_method'>
>>>
print в Python2 это statement, что можно перевести на технический русский как инструкция. print() по умолчанию т.н. “синтаксический сахар”. Т.о. это не объкт. поэтому выяснить тип не удасться -ошибка синтаксиса.
Если из __future__ модуля импортировать функцию print(), копирующая подобное поведение в Python3, то она заменит инструкцию встроенной функцией.



_________________________________________________________________________________
полезный блог о python john16blog.blogspot.com

Офлайн

#9 Окт. 7, 2015 19:47:01

Varmint
Зарегистрирован: 2015-07-23
Сообщения: 10
Репутация: +  -1  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

doza_and
ну на питоне факторизацию нет смысла писать.
Почему?
JOHN_16
print в Python2 это statement, что можно перевести на технический русский как инструкция. print() по умолчанию т.н. “синтаксический сахар”. Т.о. это не объкт. поэтому выяснить тип не удасться -ошибка синтаксиса.
Если из __future__ модуля импортировать функцию print(), копирующая подобное поведение в Python3, то она заменит инструкцию встроенной функцией.
Не знал. Спасибо!

Я совсем не сидел за компьютером пару последних дней, поэтому как только возьмусь за скрипт, так сразу напишу, что получилось, а что нет.

Офлайн

#10 Окт. 7, 2015 22:29:50

Varmint
Зарегистрирован: 2015-07-23
Сообщения: 10
Репутация: +  -1  -
Профиль   Отправить e-mail  

Решение системы уравнений методом Гаусса

Слизал с примера на Java из Википедии (пытался по алгоритму, но так и не понял, что за переменная c в индексах). Получилось вот так, но считает неправильно. Реквестирую помощь.

matrix = [[1,1,1,6],
          [1,-1,2,5],
          [2,-1,-1,-3]]
def guass(a):
	x = [a[i][-1] for i in range(len(a))]
	for k in range(1,len(a)-1):
		for j in range(k,len(a)-1):
			m = a[j][k-1]/a[k-1][k-1]
			for i in range(len(a)-1):
				a[j][i] -= m*a[k-1][i]
			x[j] -= m*x[k-1]
	for i in range(len(a)-1,-1,-1):
		for j in range(i+1,len(a)-1):
			x[i] -= a[i][j]*x[j]
			x[i] = x[i]/a[i][i]
	print(x)

Не работает return (понятия не имею, почему). Результатом является список x, который в данном примере должен состоять из элементов 1,2 и 3.

Отредактировано Varmint (Окт. 8, 2015 13:13:36)

Офлайн

Board footer

Модераторировать

Powered by DjangoBB

Lo-Fi Version