Несколько квадратичных форм

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 09 фев 2019, 23:47

Теперь о том, какой метод решения привел к успеху. Вместо исходной системы

$$ \frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}=z(T_a,\rho_{m0a};x,H) $$

рассмотрим систему с дополнительным параметром $\lambda$

$$ \frac{p_{0a}+\lambda p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}=z(T_a,\rho_{m0a};x(\lambda),H(\lambda)). $$

При $\lambda=1$ возвращаемся к исходной задаче, а при $\lambda=0$ имеем задачу

$$ \frac{p_{0a}}{10^{-3}\rho_{m0a}RT_a}=z(T_a,\rho_{m0a};x(0),H(0)), $$

решением которой является эталонная смесь: $x(0)=x_0$, $H(0)=H_0$.

Продифференцируем $\lambda$-задачу по $\lambda$

$$ \frac{p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}=\sum_{i=1}^4M_{ai}({\bf x})\frac{dx_i}{d\lambda} $$

($M_{ai}({\bf x})=\partial z(T_a,\rho_{m0a};x,H)/\partial x_i$, ради компактности обозначений я объединил $x$ и $H$ в один вектор параметров $\bf x$). Полученное дифференциальное уравнение с начальными условиями $x(0)=x_0$, $H(0)=H_0$ можно решать каким-либо известным методом, например, Эйлера или Рунге---Кутта. Значения $x(1)$, $H(1)$ составят решение исходной задачи.

Так вот, решение системы дифференциальных уравнений методом Рунге---Кутта 4-го порядка с удвоением числа шагов до достижения требуемой точности прекрасно сработало на 50000 случайно выбранных смесях (при этом даже и шагов-то интегрирования много не понадобилось, максимум 32, а в среднем, думаю, 4 или 8).

Но я, убейте меня, не понимаю, почему совершенно не работает Ньютон, у которого итерация $x\to x'$

$$ \frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}-z(T_a,\rho_{m0a};x,H)= \sum_{i=1}^4 M_{ai}({\bf x})(x_i'-x_i) $$

выполняется с той же самой матрицей $M$! Причем я пытался скрестить два метода: делал маленькие шаги по $\lambda$ и каждый раз пытался дотянуть решение Ньютоном с предыдущего шага --- так фиг тебе! Не дотягивает ни при каком разумном шаге.

Вот причину этих странностей мне и хотелось бы обсудить. Рабочая гипотеза у меня есть: дело в том, что нелинейность мала, а линейная часть уравнений вырождена. Вырождение это не случайное, о нем говорилось еще в заглавном сообщении темы: для смеси идеальных газов восстановить состав по измерениям $p,V,T$ невозможно.

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 10 фев 2019, 00:55

peregoudov писал(а):Source of the post а линейная часть уравнений вырождена

Это легко проверить.
Проанализируйте численно полученную матрицу на обусловленность - там видно будет.

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 11 фев 2019, 10:00

Вырожденность линейной части и проверять нечего: от $H$ она вообще не зависит, а от $x$ зависит только в комбинации $x_e+x_a+x_y$. Так что 4х4 матрица имеет ранг 1. Одно ее собственное значение (соответствующее собственному вектору (1,1,1,1)) равно 3, остальные равны нулю.

А полная матрица... Ну вот, например, тот же случай $(T,\rho_{m0})=(290,10/22.4)$, (310, 10/22.4), (312, 1/22.4), (330, 1/22.4), $x_{\textrm{e}}=0.93$, $x_{\textrm{a}}=0.05$, $x_{\textrm{y}}=0.02$, $\rho_c=0.86$ ($H=1065.33$), что на предыдущей странице. Сама матрица производных имеет вид

$$ \begin{pmatrix} 0.937790&0.975811&0.926953&-7.068753\cdot10^{-5}\\ 0.947311&0.981041&0.937826&-6.141893\cdot10^{-5}\\ 0.994626&0.998013&0.993683&-0.619812\cdot10^{-5}\\ 0.995354&0.998368&0.994480&-0.548691\cdot10^{-5} \end{pmatrix}. $$

А ее собственные значения

2.897145, 0.0152640, 0.992192e-4, 0.397268e-6.

Насколько я понимаю, число обусловленности --- это просто отношение максимального и минимального собственных значений. Ну так оно порядка $10^7$.

Аватар пользователя
Ian
Сообщений: 960
Зарегистрирован: 18 янв 2016, 19:42

Несколько квадратичных форм

Сообщение Ian » 11 фев 2019, 18:11

Метод Ньютона это не метод"пытки с помощью веревочной петли и палки", в котором, как только сказано название, сразу понятно, что делают. Пусть [math]-вектор, минимизируется [math]. Определяется направление смещения [math], где G матрица Гессе f в точке x. Модифицированный метод Ньютона состоит в том, что смещаемся на [math], где [math] определяется из одномерной задачи[math]. Далее, G может оказаться необратимой (вырожденной), или операция обращения связана с большими численными ошибками. Поэтому придумана группа "Квазиньютоновских методов", которые избегают обращения матрицы G на каждом шаге итераций. Действительно, [math] непрерывно зависит от х,и изменение матрицы [math], при малом (мы надеемся) d можно линеаризовать, и тоже довольно правильно будет давать направление итераций. Я расскажу несколько методов, если надо. У Вас оптимизация по трем переменным?

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 12 фев 2019, 09:39

Ian, мне снова кажется, что вам будет значительно легче понимать, о чем идет речь, если вы попытаетесь прочитать, что пишут другие.

Во-первых, мы не минимизаций функции занимаемся, а решением системы уравнений.

Во-вторых, неизвестных не три, а четыре (и уравнений тоже четыре).

В-третьих, если вы прочтете первое сообщение на этой странице, вы увидите, что там приведены итерационные шаги для Ньютона, который работать начисто отказывается, и для $\lambda$-метода, который прекрасно работает. Оба этих метода используют одну и ту же матрицу ($M$ в моих обозначениях и $G$ в ваших, кстати, та же матрица возникает в задаче оценки погрешности решения уравнения по погрешностям измерений), слегка различаясь лишь правыми (в моих формулах --- левыми) частями. Вот эта разительная разница в поведении и есть самая большая непонятка.

Аватар пользователя
Ian
Сообщений: 960
Зарегистрирован: 18 янв 2016, 19:42

Несколько квадратичных форм

Сообщение Ian » 12 фев 2019, 17:10

peregoudov писал(а): попытаетесь прочитать, что пишут другие.
Именно что пытаюсь.Ну меня -то вы поняли, судя по ответу. Так кто тут молодец, кто пишет понятно или кто непонятно :)

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 13 фев 2019, 09:57

Давайте сойдемся на том, что молодец --- тот кто понимает даже непонятно написанное :) Вы своими мыслями поделитесь, отчего так странно получается? И как преодолеть плохую определенность матрицы?

Аватар пользователя
Ian
Сообщений: 960
Зарегистрирован: 18 янв 2016, 19:42

Несколько квадратичных форм

Сообщение Ian » 13 фев 2019, 13:56

Вопрос с 99,99% уже изучен всесторонне только надо поставить его как частный случай известной задачи. Я знаю все известные задачи, но не знаю какая из них у вас. Вы наоборот не знаете известных задач почти ни одной(

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 13 фев 2019, 15:47

Так для того и поднята опять тема, чтобы разобраться, что у меня за задача. С реальными задачами всегда ведь так: самое трудное --- понять, какая из стандартных задач тут возникает. Ну, и что вам нужно, чтобы это понять? Вот zykov сперва просил в явном виде написать уравнения. Ну, я написал, вам из них понятнее стало? Потом он про матрицу спросил --- вот вам матрица, вот вам собственные значения, какие выводы вы можете сделать?

Могу я попытаться угадать, какой стандартной задаче это все соответствует. Мне представляется что-то типа

$$ b_i+\varepsilon c_i+\sum_j A_{ij}x_j+\varepsilon f_i({\bf x})=0,\quad i,j=1,\ldots,n, $$

причем матрица $A$ вырождена, но уравнение $Ax=b$ имеет решение(я). Тогда при $\varepsilon=0$ как раз имеем это решение, а при малых $\varepsilon$...

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 13 фев 2019, 19:50

peregoudov писал(а):Source of the post причем матрица $A$ вырождена, но уравнение $Ax=b$ имеет решение(я)

Я с таким не сталкивался, но уверен, что люди давно придумали методы.
Наверно нужно привести задачу (переформулировать) к такому виду, чтобы матрица была не вырожденной (например уменьшая размерность).

peregoudov писал(а):Source of the post Ну так оно порядка $10^7$.

Не знаю, насколько Вы знакомы с вычислительными методами (ошибки округления машинных чисел и прочее).
Для числа с плавующей точкой одинарной точности (single float) $10^7$ - это много. Там как раз относительная ошибка округления порядка $10^{-7}$.
Для двойной точности (double float) относительная ошибка округления порядка $10^{-16}$, так что должно быть достаточно. Но если например где-то возникает квадрат $10^{-7}$, то уже могут быть проблемы.

Один из простейших подходов, который используют при вычислениях - это нормировка. Работают не с начальными значениями, а с нормированными. Например если одни из величин имеют порядок $10^2$ (например в вольтах), а другие порядок $10^{-5}$ (например в амперах), то и те, и другие нормируют, чтобы они имели порядок $10^0$.

Вообще, составление вычислительной схемы - это отдельная большая работа (нужно следить за ошибками, за стабильностью и прочее). Если у вас нет таких знаний, то нужно либо серьёзно это изучать, либо привлекать другого человека, кто разбирается в этих вопросах.

Аватар пользователя
Ian
Сообщений: 960
Зарегистрирован: 18 янв 2016, 19:42

Несколько квадратичных форм

Сообщение Ian » 13 фев 2019, 20:07

я пытаюсь понять но все равно все запутано. стереть бы тему и начать сначала более мелкими шагами

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 13 фев 2019, 20:14

zykov писал(а):Source of the post либо привлекать другого человека, кто разбирается в этих вопросах
Ну, а я чем сейчас занимаюсь?

Масштабирование пробовал, не работает. А вы уверены, что оно должно работать для Ньютона, а не просто тупо сокращаться?

Считается все по умолчанию в С, с двойной точностью.

Я далеко не обо всем написал, что я перепробовал за эти полгода... Например, я теперь умею пользоваться gmp.

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 13 фев 2019, 20:15

Ian писал(а):Source of the post стереть бы тему и начать сначала более мелкими шагами
Мелкими шагами --- это можно, а стирать-то зачем?

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 13 фев 2019, 23:33

peregoudov писал(а):Source of the post Ну, а я чем сейчас занимаюсь?

Имелось ввиду, что человек нанимается за соответствующую плату для выполнения определенного объема работ.

peregoudov писал(а):Source of the post Считается все по умолчанию в С, с двойной точностью

По умолчанию в C тип float имеет одинарную точность. Двойную точность имеет тип double.

peregoudov писал(а):Source of the post за эти полгода

Полгода это срок...

А можно на пальцах объяснить, откуда берется вырожденность матрицы? Это именно принципиальная вырожденность или просто плохообусловленность?
Просто если матрица вырождена, то это не лучший способ описать требуемое решение. Нужно как-то по другому описать то, что нужно найти...

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 14 фев 2019, 00:42

peregoudov писал(а):Source of the post Но я, убейте меня, не понимаю, почему совершенно не работает Ньютон, у которого итерация $x\to x'$

$$ \frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}-z(T_a,\rho_{m0a};x,H)= \sum_{i=1}^4 M_{ai}({\bf x})(x_i'-x_i) $$

Здесь случайно знак не перепутан? А то может шаг делается в противоположенном направлении...
Например в одномерном случае шаг делается исходя из $f(x)=f'(x)(x-x')$ (тут $x'$ - новое значение).

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 14 фев 2019, 19:24

zykov писал(а):Source of the post Имелось ввиду, что человек нанимается за соответствующую плату для выполнения определенного объема работ.
Работы выполнены, заказчик удовлетворен. Вы, видимо, не совсем понимаете моих побудительных мотивов. Дело не только в том, что они мне платят деньги (и даже не столько в том: если посчитать интегрально, сколько они заплатили за год и поделить помесячно --- выйдет очень смешная сумма). Я во всем пытаюсь найти какой-то академический интерес, посмотреть на задачу шире, чем она изначально поставлена, узнать в процессе решения что-то новое. Просто сговнякать что-то непонятно почему хоть как-то работающее --- мне не интересно (хотя вполне удовлетворяет заказчика).

zykov писал(а):Source of the post По умолчанию в C тип float имеет одинарную точность. Двойную точность имеет тип double.
Основным типом чисел с плавающей точкой в С является именно double. Например, sin(x) --- это функция, принимающая и возвращающая тип double. Соответствующая функция для типа float называется sinf(x).

zykov писал(а):Source of the post А можно на пальцах объяснить, откуда берется вырожденность матрицы? Это именно принципиальная вырожденность или просто плохообусловленность?
Я вроде объяснял уже несколько раз... Дело в том, что в уравнение идеального газа не входит никаких индивидуальных характеристик этого газа (вообще у идеального газа всего одна такая характеристика --- теплоемкость, и она входит во внутреннюю энергию). Поэтому состав смеси идеальных газов невозможно определить по измерениям p, V, T. Тут вырожденность точная.

Уравнение GERG91, которое я привел --- это смесь слабонеидеальных газов, то есть малые поправки к уравнению идеального газа. Тут матрица получается невырожденной, но плохо обусловленной (тем хуже, чем ближе смесь к идеальному газу).

zykov писал(а):Source of the post Здесь случайно знак не перепутан?
Да вроде нет. Линеаризуется уравнение

$$ \frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}=z(T_a,\rho_{m0a};x',H') $$

вблизи $(x,H)$, что получается в правой части? Вроде как

$$ z(T_a,\rho_{m0a};x,H)+\sum_{i=1}^4 M_{ai}({\bf x})(x_i'-x_i). $$

Первое слагаемое переносим в левую.

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 14 фев 2019, 21:54

peregoudov писал(а):Source of the post Да вроде нет

Если матрица $M_{ai}$ соответствует производной, то точно перепутан. Если минус-производной, то всё верно.
(Нигде не нашел определения этой матрицы.)
Впрочем в программе проверить легко. Можно поменять знак и посмотреть, начнет ли сходится.

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 14 фев 2019, 22:39

peregoudov писал(а):Source of the post Я во всем пытаюсь найти какой-то академический интерес

Не в обиду сказано.
Со стороны выглядит, как "вот несколько страниц кривого кода, помогите его отладить".

peregoudov писал(а):Source of the post то есть малые поправки к уравнению идеального газа. Тут матрица получается невырожденной, но плохо обусловленной

Если матрица имеет вид "вырожденная + малая поправка", то можно попробовать найти способ решения для этого случая с контролем вычислительной ошибки.

Проблема с плохообусловленной матрицей в том, что малая ошибка в правой части (а она там всегда есть в силу округления машинного числа) приводит к большой ошибке в решении.
Можно например попробовать так. Если $x_0$ - решение невозмущенной системы $A_0 x_0=b_0$, то искать решение возмущенной системы как $(A_0+A_1)(x_0+x_1)=b_0+b_1$, где $A_1$, $x_1$, $b_1$ малы. Важно, что малые величины вычисляются напрямую с полной точностью, а не как разность (например $b_1$ вычисляется напрямую, а не как разность двух больших $(b_0+b_1)-b_0$). Тогда $(A_0+A_1)x_1=b_1-A_1 x_0$. В правой части оба слагаемых малы и имеют полную точность. Найденное $x_1$ конечно не будет иметь полной точности в силу плохой обусловленности матрицы, но будет более точным, чем при вычислении с большой правой частью.
При необходимости можно сделать несколько итераций в этом направлении.

peregoudov
Сообщений: 620
Зарегистрирован: 29 дек 2015, 13:17

Несколько квадратичных форм

Сообщение peregoudov » 15 фев 2019, 12:11

zykov писал(а):Source of the post (Нигде не нашел определения этой матрицы.)
Оно здесь
peregoudov писал(а):Source of the post$M_{ai}({\bf x})=\partial z(T_a,\rho_{m0a};x,H)/\partial x_i$,
Никакой ошибки в знаке по-прежнему не вижу. Потом, в программе же используется готовая процедура решения Ньютоном, туда загоняются процедура вычисления значения самой функции и ее производной, важно, чтобы они были согласованы по знаку. Проверил еще раз --- согласованы.

zykov писал(а):Source of the post Со стороны выглядит, как "вот несколько страниц кривого кода, помогите его отладить".
Тот код, что я прикреплял летом, он же состоит из частей. gerg91plus.c --- это просто запрограммированные формулы из ГОСТа. Эта часть проверялась параллельным программированием в Maple и сравнением результатов тестовых расчетов, кроме того, сравнивалась с другой реализацией метода gerg91 (код дали мои нефтяники). Ну то есть эту часть отлаживать не требуется.

Кучка файлов fdjac.c, fmin.c, lnsrch.c, lnbksb.c, ludcmp.c, newt.c, nrutil.c --- это все готовый код из Numerical Recipes, его тоже отлаживать не надо, хотя я всегда проверяю на паре тестов перед первым использованием.

Остается один файл solve_gerg.c, в котором всего несколько строк, в частности, именно там заданы процедуры вычисления левой части уравнений и производных от нее. Напортачить можно было только здесь.

zykov писал(а):Source of the post Можно например попробовать так. Если $x_0$ - решение невозмущенной системы $A_0 x_0=b_0$, то искать решение возмущенной системы как $(A_0+A_1)(x_0+x_1)=b_0+b_1$, где $A_1$, $x_1$, $b_1$ малы.
Идея понятная. Вообще-то полное уравнение нелинейное, но главное не в этом. Вы исходите из ошибочного предположения, что $x_0$ единственно. Нет, тут как в квантовой механике при вырождении уровня: в нулевом порядке вы можете лишь сказать, что $x_0$ лежит в некотором --- неодномерном! в нашем случае размерности 3 --- подпространстве. Конкретное значение определится лишь в следующем приближении. Я эту теорию возмущений крутил и даже программировал, но точность первого (или первых двух, не помню уже точно) приближения недостаточна.

zykov
Сообщений: 1393
Зарегистрирован: 06 янв 2016, 17:41

Несколько квадратичных форм

Сообщение zykov » 16 фев 2019, 01:12

peregoudov писал(а):Source of the post $$ \frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}-z(T_a,\rho_{m0a};x,H)= \sum_{i=1}^4 M_{ai}({\bf x})(x_i'-x_i) $$
Судя по этому матрица $M_{ai}$ должна соответствовать минус производной $\frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}-z(T_a,\rho_{m0a};x,H)$. Если производная от $\frac{p_{0a}+p_{\Delta a}}{10^{-3}\rho_{m0a}RT_a}$ равна нулю, то так оно и есть.

Я бы советовал занятся отладкой, тогда может что-то прояснится.
Во-первых отказатся от сторонней процедуры для Ньютона (кто знает, что они туда ещё заложили). А написать этот алгоритм самостоятельно (это элементарно, когда уже есть и значение функции, и значение производной).
Во-вторых прогнать несколько итераций этого алгоритма и записать в лог-файл поведение алгоритма. Какие $x$ на каждой итерации и какой квадрат ошибки целевой функции.
Если точное решение известно, то будет видно как себя ведет $x$ по итерациям - приближается ли, удаляется ли, стоит на месте или вообще скачет. И по квадрату ошибки целевой функции будет видно, уменьшается ли он, и если да, то как быстро.

Вообще метод Ньютона - это градиентный спуск для квадрата ошибки. Плюс ещё, Ньютон даёт не только направление движения, но и размер шага. Но размер шага хорош только когда нелинейность невелика. Поэтому в этом случае делают маленькие шаги, т.е. делают просто спуск по граденту.
Если посмотреть на трейс этого квадрата ошибки то будет видно, работает алгоритм или нет. Если ошибок нет, то этот квадрат должен на каждом шаге убывать (т.к. спуск по градиенту) и в пределе стремится к нулю.


Вернуться в «Математика»

Кто сейчас на форуме

Количество пользователей, которые сейчас просматривают этот форум: нет зарегистрированных пользователей и 12 гостей