DjangoBB LoFi version

Полная версия: как написать свой солвер для систем нелинейных уравнений?

Начало » Python для экспертов » как написать свой солвер для систем нелинейных уравнений?

somesanctus

Июнь 8, 2011 10:34:32

на настоящий момент пользуюсь функционалом nsolve из SymPy, но не хочется ради одного уравнения использовать дополнительные библиотеки.
система трех уравнений остается неизменной, меняются только 9 коэффицентов, первое приближение тоже всегда одинаковое

from sympy import *

A=[-297.931,0,-56.708]
B=[-30.389,0,1263.047]
C=[219.725,0,-42.816]

def findPosition(A,B,C):
	x,y,z = symbols('xyz')

	x1,y1,z1=A
	x2,y2,z2=B
	x3,y3,z3=C

	# equations
	S1=(x-(x1+x2)/2)**2+(y-(y1+y2)/2)**2+(z-(z1+z2)/2)**2-(sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)/2)**2
	S2=(x-(x2+x3)/2)**2+(y-(y2+y3)/2)**2+(z-(z2+z3)/2)**2-(sqrt((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)/2)**2
	S3=(x-(x3+x1)/2)**2+(y-(y3+y1)/2)**2+(z-(z3+z1)/2)**2-(sqrt((x3-x1)**2+(y3-y1)**2+(z3-z1)**2)/2)**2

	answer = nsolve((S1,S2,S3),(x,y,z),(0,1,0))
	return answer
print findPosition(A,B,C)

как написать свой солвер, с чего начать, где подсмотреть реализации?

doza_and

Июнь 8, 2011 20:26:06

:) см раздел изобретение велосипедов. На этом форуме нет, наверно надо завести.
http://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9D%D1%8C%D1%8E%D1%82%D0%BE%D0%BD%D0%B0
решалки нелинейных уравнений можно найти еще в numpy и посмотреть там исходный код.

somesanctus

Июнь 9, 2011 02:16:16

есть ли здесь математики?)))

по алгоритму ньютона в принципе понятно

берем первое приближение, находим значение функции в нем, делим на значение производной в этом же приближении, получаем следующее приближение, повторяем пока не покроем погрешность.

но не понятно как работать с системами уравнений. что-то подсказывает что они должны стать элементами матрицы и все решаться будет именно ими.
прав ли я?

doza_and

Июнь 9, 2011 10:08:34

somesanctus
должны стать элементами матрицы

Да вы правы это один из множества способов. И все это уже давным давно написано. Писать свою библиотеку мотивируя это тем что вам не хочется тащить с собой библиотеку - крайне странно. У вас будет все равно библиотека, только своя :). Численное решение потащит за собой матричную алгебру и оценку машинной точности, Алгоритмы решения некорректных задачек(для дифференцирования) алгоритмы локализации корней, систему сообщения об ошибках и может еще чего. Если нужна производительность - надо будет освоить изготовление динамических библиотек и язык которым бы их будете делать кроме того надо будет делать систему автоматического дифференцирования - для аналитического или полуаналитического расчета матрицы Якоби. Свой словер я могу понять толлько если вы придумали хороший прием для решения именно ваших уравнений и вам требуется очень быстрый решатель.
посмотрите простенькие книжки Калиткин Н.Н.: “Численные методы”, Москва: Наука, 1978
Самарский А.А.: “Введение в численные методы.”, Москва: Наука, 1978.
Фаддеев Д.К., Фаддеева В.Н.: “Вычислительные методы линейной алгебры.”, Москва: ФМЛ, 1960.
Мину М.: “Математическое программирование. Теория и алгоритмы.”, М.: Наука., ФМЛ, 1990
Для более подробного ознакомления с вопросом.
Трауб Дж.: “Итерационные методы решения уравнений”, М.: МИР, 1985

Если грубо то так:
f'(x)^(-1)->matrix(df_i/dx_j)^(-1)

somesanctus

Июнь 9, 2011 10:35:41

спасибо за подробности.
собственно я написал плагин для MAYA, который “восстанавливает камеру”
посмотреть его в действии можно здесь http://vimeo.com/24130335
CG сообщество заинтересовалось, но когда я написал, что потребуется установка дополнительных либ, желающих потестить заметно поубавилось. тем более что под виндой все не так облачно с этим, как в nix и первым пожеланием коллег была самодостаточность плагина.

собственно почему я решил, что задача не так уж и сложна:
1. уравнение неизменно, т.е.:
2. пересчитывать для него производные не нужно,
3. первое приближение мне всегда известно, хотябы визуально, на деле же небольшого отступа от нуля хватает, чтобы получать значение с двойной точностью за 0,02 сек(вполне удовлетворяет работу в реалтайме).

мне кажется, что учитывая это, большую часть того, что вы написали в контексте “потянет за собой” можно опустить.

ZAN

Июнь 9, 2011 12:35:34

Ничего сложного там нет, если о конкретике, то ссылка на конспект здесь:
http://physlib.comule.com/book/Proga_Shapochkina_konsp_3sem.zip

А здесь уже без каракулей
http://www.nsu.ru/matlab/Exponenta_RU/educat/systemat/hanova/equation/loc.asp.htm

doza_and

Июнь 9, 2011 13:00:37

Да вы правы, можно все упростить, но помоему правильный путь всетаки сделать дистрибутивчик и не беспокоится. Вы получите решение которое работает в 99% случаев а потом можете 3 года ловить частные случаи. Если задачка хорошо сходится то не парьтесь в ньютоном сделайте простые итерации и все. Если плохо сходятся то можно нарваться на потерю сходимости в ньютоне.