Форум сайта python.su
Пытаюсь написать скрипт, который бы решал методом Гаусса матрицы хотя бы до 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)
Отредактировано Varmint (Окт. 6, 2015 20:15:29)
Офлайн
VarmintМда. А что искали?
Искал в Интернете готовые скрипты - не подходят даже после каких-либо изменений.
Офлайн
VarmintВообще, через __future__ во втором питоне можно импортировать функцию print().
P.S. Пусть вас не смущает функция print. Я пишу на Pythonista. Там разработчик предусмотрел возможность писать print со скобочками и без.
>>> from __future__ import print_function >>> >>> print(1, 2, 3, sep='x') 1x2x3 >>>
VarmintДавным-давно я делал для NxM, но так и не сделал. Сделал только приведение к треугольному виду, но система может иметь фундаментальную систему решений (когда для точных значений переменных не хватает информации).
Пытаюсь написать скрипт, который бы решал методом Гаусса матрицы хотя бы до 10х10 включительно.
Отредактировано py.user.next (Окт. 6, 2015 02:48:08)
Офлайн
doza_andСуть как раз-таки в том, что надо самому написать, а с NumPy я не сталкивался ещё, поэтому даже сурс плохо понимаю. Если Вы можете переписать использованные в примере функции на обычный язык, то огромное Вам спасибо.
Мда. А что искали?
http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html
py.user.nextВ моём случае это только квадратная.
А смысл делать алгоритм для квадратной, если основная матрица системы может быть неквадратной?
Офлайн
# 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
Офлайн
ну на питоне факторизацию нет смысла писать.
Простейший код выглядит примерно так:
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]; } } } }
#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
Отредактировано doza_and (Окт. 6, 2015 21:42:46)
Офлайн
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)
Офлайн
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'> >>>
Офлайн
doza_andПочему?
ну на питоне факторизацию нет смысла писать.
JOHN_16Не знал. Спасибо!
print в Python2 это statement, что можно перевести на технический русский как инструкция. print() по умолчанию т.н. “синтаксический сахар”. Т.о. это не объкт. поэтому выяснить тип не удасться -ошибка синтаксиса.
Если из __future__ модуля импортировать функцию print(), копирующая подобное поведение в Python3, то она заменит инструкцию встроенной функцией.
Офлайн
Слизал с примера на 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)
Отредактировано Varmint (Окт. 8, 2015 13:13:36)
Офлайн