Модифицированный метод Грахама для нахождения начальных значении для решения системы нелинейных уравнении методом Ньютона-Рафсона из книги
Наташка Перевощикова
1)Три вектора:
a) вектор с site fraction по углероду в интервале от 9.е-7 до
9.е-1.Можно провести дисцритизацию по мольной доли (Xc) в интервале от
0 до 1, если вы работаете в уравнениях с мольнями долями. Если же в
уравнениях для энергии Гиббс работаете с site fraction, то
дисцритизацию надо проводить по Yc.При этом в интервале не может
фигурировать от 0 до 1 , так как Gmix = 1/a‘yi ln(yi)
б) вектор с энергией Гиббс по фазе номер 1 в зависимости от Yc
г) вектор с энергией Гиббс по фазе 2 в зависимости от
2)Сравниваем два вектора с энергией Гиббс по фазе 1 и пофазе 2, и
выбираем минимальные значения. Таким образом получаем область нижней
части энергии Гиббс - область выпуклой оболочки:
G[n]<-min(Gph1,Gph2)
3) if (G[n]==Gph1)
{
iph[n]=0; //Фаза номер 1
}
else
{
iph[n]=1 //Фаза номер 2
}
4) Graham scan modifie
Теперь выбираем точки для выпуклой оболочки
ihull[n];
double ihullmax = 0.0;
for (i = 2; i < n; i++)
{
.В качестве условия построения выпуклой оболочки использую это
(Это для двух компонентной системы железо углерод):
∆Gf(P,T,x)<= [(x“-x)/(x”-x‘)]∆Gf(P,T,x‘)+[(x-x‘)/(x”-x‘)]∆Gf
(P,T,x”)
function ∆Gf(Y[ihullmax], Y[i-1], Y[i]);
if ( function ∆Gf >= 0.0)
{
ihullmax++;
ihull[ihullmax]=i-1;
}
else
{
не изменяется ihullmax и берется следущie tochki для
проверки на функцию ∆Gf...Таким образом должны выбрасываться не нужные
точки, если конечно я использую верное условие.
}
ihull.resize(ihullmax);
for (i = 0; i < ihullmax; i++)
{
if (iph[ihull[i]]!=iph[ihull[i-1]]) return
Yc для фазы 1 = Y[ihull[i-1]] ;
Yc для фазы 2 = Y[ihull[i]];
}
Вот и весь алгоритм по нахождению начальных значении двух компонентной
системы железо-углерод. Должна возврашать в качестве значении site
fraction углерода, где имеется гетерогенное равновесное состояние
Теперь:
Возвращает значения:
Yc phase1 =2.698900e-03 ;
Yc phase2 =3.046900e-03 ;
В интервале ot 9.e-7 do 4.e-2 с шагом 1.e-6:
convexhull(9.e-7,4.e-2,1.e-6,y1alfa,y2alfa,y1betta,y2betta);
Thermocalc:
W(c%)=2.e-3;
Y(BCC_A2,C#2)=2.9488045E-4
Y(FCC_A1,C#2)=3.558066E-2
Это не близкие значения к термокальки и поэтому Ньютон-Рафсон не
работает без медленного спуска
x[i] = xold[i] + alam * p[i]* 1e-3;
Не знаю, где ошибка!
Какой именно метод Ньютона лучше всего использовать в этом случае?
Evgenii Rudnyi
Давайте последовательно. Вначале общие соображения.
> Не знаю, где ошибка!
Ошибка может быть везде. Здесь надо вначале построить мольные энергии
Гиббса как функции от мольной доли или site fraction. Тогда можно будет
визуально увидеть правильный ответ. Если он не совпадает с результатами
вашей программы, то надо делать промежуточный вывод, шаг за шагом, и
сравнивать с тем, что должно быть. Вполне возможно, что алгоритм
правильный, но где-нибудь есть какая-то глупая ошибка при
программировании (скорее всего и не одна).
Теперь более конкретные вопросы.
Итак, задача это расчет равновесного состава в системе Fe-C. Какие фазы
принимаются в расчет? Какие модели фаз используются? Вы берете энергии
Гиббса из Thermocalc, или они запрограммированы вами?
> 1)Три вектора:
> a) вектор с site fraction по углероду в интервале от 9.е-7 до
> 9.е-1.Можно провести дисцритизацию по мольной доли (Xc) в интервале от
> 0 до 1, если вы работаете в уравнениях с мольнями долями. Если же в
> уравнениях для энергии Гиббс работаете с site fraction, то
> дисцритизацию надо проводить по Yc.При этом в интервале не может
> фигурировать от 0 до 1 , так как Gmix = 1/a‘yi ln(yi)
> б) вектор с энергией Гиббс по фазе номер 1 в зависимости от Yc
> г) вектор с энергией Гиббс по фазе 2 в зависимости от
> 2)Сравниваем два вектора с энергией Гиббс по фазе 1 и пофазе 2, и
> выбираем минимальные значения. Таким образом получаем область нижней
> части энергии Гиббс - область выпуклой оболочки:
> G[n]<-min(Gph1,Gph2)
> 3) if (G[n]==Gph1)
> {
> iph[n]=0; //Фаза номер 1
> }
> else
> {
> iph[n]=1 //Фаза номер 2
> }
То есть из G[n] вы теперь хотите построить выпуклую оболочку?
> 4) Graham scan modifie
> Теперь выбираем точки для выпуклой оболочки
> ihull[n];
> double ihullmax = 0.0;
Этот кусочек пока остался непонятным. А почему вы не используется уже
готовые подпрограммы для нахождения выпуклой оболочки?
> Вот и весь алгоритм по нахождению начальных значении двух компонентной
> системы железо-углерод. Должна возврашать в качестве значении site
> fraction углерода, где имеется гетерогенное равновесное состояние
> Теперь:
> Возвращает значения:
>
> Yc phase1 =2.698900e-03 ;
> Yc phase2 =3.046900e-03 ;
> В интервале ot 9.e-7 do 4.e-2 с шагом 1.e-6:
> convexhull(9.e-7,4.e-2,1.e-6,y1alfa,y2alfa,y1betta,y2betta);
>
> Thermocalc:
> W(c%)=2.e-3;
> Y(BCC_A2,C#2)=2.9488045E-4
> Y(FCC_A1,C#2)=3.558066E-2
Как я написал выше, надо отлаживать программу. Возможны ошибки
программирования, возможно дискретизация слишком грубая, возможны
логические ошибки. Грубо говоря, вам надо прогнать ваш алгоритм один раз
вручную, то есть, делать вывод и проверять его после каждого шага.
> Это не близкие значения к термокальки и поэтому Ньютон-Рафсон не
> работает без медленного спуска
> x[i] = xold[i] + alam * p[i]* 1e-3;
> Не знаю, где ошибка!
> Какой именно метод Ньютона лучше всего использовать в этом случае?
Я думаю, что универсального метода нет. Я в свое время использовал
DONLP2, поскольку эта подпрограмма учитывает ограничения. Но хорошее
начальное приближение там также требовалось.
Наташка Перевощикова
Гиббса как функции от мольной доли или site fraction. Тогда можно будет
визуально увидеть правильный ответ. Если он не совпадает с результатами
вашей программы, то надо делать промежуточный вывод, шаг за шагом, и
принимаются в расчет? Какие модели фаз используются? Вы берете энергии
Гиббса из Thermocalc, или они запрограммированы вами?
готовые подпрограммы для нахождения выпуклой оболочки?
DONLP2, поскольку эта подпрограмма учитывает ограничения. Но хорошее
начальное приближение там также требовалось.
Evgenii Rudnyi
> python.Результаты
> проверки сходятся.Затем построила с помошью опять таки python кривые
> зависимости энергии Гиббс от site fraction углерода.Визуально увидела
> правельный ответ.Модель энергий гиббс должна отражать модель thermocalc, так
> как я делаю копию thermocalc и сравниваю каждый раз свой результаты с
> thermocalc. Результатами сравнения я осталась довольна
Это уже непонятно. Какие результаты показывает графическое постороение?
Совпадают ли они с Thermocalc? Или нет?
То есть, что означает "Результатами сравнения я осталась довольна"? Есть
разница с Thermocalc, или нет?
> >Итак, задача это расчет равновесного состава в системе Fe-C. Какие
> фазы
> принимаются в расчет? Какие модели фаз используются? Вы берете
> энергии
> Гиббса из Thermocalc, или они запрограммированы вами?
> Расчет равновесного состава в системе железо-углерод (феррит -
> цементит).Энергии Гиббс запрагроммировала сама и сравнила с thermocalc.
Модель фазы - регулярный раствор?
> >То есть из G[n] вы теперь хотите построить выпуклую оболочку?
> Этот вектор представляет собой в геометрическом плане нижнюю часть профиля
> энергии гиббс, которая нас очень интересует. Из этого вектора нужно выбрать
> корретктно точки для построения выпуклой поверхности
То есть, в вашем случае две точки, которые надо соединить прямой линией,
которая будет описывать гетерогенную область.
> >Этот кусочек пока остался непонятным. А почему вы не используется уже
> готовые подпрограммы для нахождения выпуклой оболочки?
> Хороший вопрос. Если бы я была уверенна на 120% процентов, что они "готовы"
> рассчитать многокомпонентные системы, то,навеное бы и использовала. Я начала
> с самой простой системы, а потом должна перейти к системе из 7 компонентов.
> Есть уверенность того, что эти программы будут работать, когда ты не знаешь,
> что у них в нутри? К тому же , мне самой интересно разобраться в этом. Я и
> программ таких не знаю.
Например
Наверное, если поискать, то можно найти еще.
Я согласен, что полезно знать, как такие программы работают. Но здесь
тогда надо разделить задачи. Написание общей программы построения
выпуклой оболочки - это своя работа. Знакомство с готовыми программами
ускорит разработку. Там по крайней мере есть тесты, на которых можно
проверить работу своей программы. Также полезно сравнить
производительность.
> >Я думаю, что универсального метода нет. Я в свое время использовал
> DONLP2, поскольку эта подпрограмма учитывает ограничения. Но
> хорошее
> начальное приближение там также требовалось.
> Не могли бы Вы дать полное название без этой сложной абривиатуры, потому что
> я такого не встречала. И в чем этот метод заклучается и есть ли у него
> преимушества перед методом Ньютона-Рафсона где используется функцуя поиска
> lnsrch(....)?
Основное преимущество, что вы задаете ограничения, которые будут
использоваться во время оптимизации. Традиционный алгоритм
Ньютона-Рафсона предполагает, что поиск ведется во всем пространстве
значений, что вообще говоря не так в случае термодинамике.
Про DONLP2 (Do nonlinear programming) можно посмотреть на сайте его автора
http://www.mathematik.tu-darmstadt.de:8080/ags/ag8/Mitglieder/spellucci_de.html
Разные алгоритмы минимизации есть в
Decision Tree for Optimization Software
http://plato.la.asu.edu/guide.html
А Ньютон-Рафсон вы также сами программировали?
Наташка Перевощикова
> Программу я пишу на С++. Аналитически проверку провожу вЭто уже непонятно. Какие результаты показывает графическое постороение?
> python.Результаты
> проверки сходятся.Затем построила с помошью опять таки python кривые
> зависимости энергии Гиббс от site fraction углерода.Визуально увидела
> правельный ответ.Модель энергий гиббс должна отражать модель thermocalc, так
> как я делаю копию thermocalc и сравниваю каждый раз свой результаты с
> thermocalc. Результатами сравнения я осталась довольна
Совпадают ли они с Thermocalc? Или нет?
То есть, что означает "Результатами сравнения я осталась довольна"? Есть
разница с Thermocalc, или нет?
> >Итак, задача это расчет равновесного состава в системе Fe-C. КакиеМодель фазы - регулярный раствор?
> фазы
> принимаются в расчет? Какие модели фаз используются? Вы берете
> энергии
> Гиббса из Thermocalc, или они запрограммированы вами?
> Расчет равновесного состава в системе железо-углерод (феррит -
> цементит).Энергии Гиббс запрагроммировала сама и сравнила с thermocalc.
> >То есть из G[n] вы теперь хотите построить выпуклую оболочку?То есть, в вашем случае две точки, которые надо соединить прямой линией,
> Этот вектор представляет собой в геометрическом плане нижнюю часть профиля
> энергии гиббс, которая нас очень интересует. Из этого вектора нужно выбрать
> корретктно точки для построения выпуклой поверхности
которая будет описывать гетерогенную область.
Именно
Да, сама
Evgenii Rudnyi
>>> python.Результаты
>>> проверки сходятся.Затем построила с помошью опять таки python кривые
>>> зависимости энергии Гиббс от site fraction углерода.Визуально увидела
>>> правельный ответ.Модель энергий гиббс должна отражать модель thermocalc,
>> так
>>> как я делаю копию thermocalc и сравниваю каждый раз свой результаты с
>>> thermocalc. Результатами сравнения я осталась довольна
>> Это уже непонятно. Какие результаты показывает графическое постороение?
>> Совпадают ли они с Thermocalc? Или нет?
>
> Графические результаты совпадают с термокалькой
>
>>
>> То есть, что означает "Результатами сравнения я осталась довольна"? Есть
>> разница с Thermocalc, или нет?
>>
> Разницы нет, одинаковые
Но тогда это означает, что ошибка в программе. То есть ваша программа
выдает значения, которые расходятся с графическим построением. Разве не
так? Извините, но сейчас стало совсем непонятно.
>> А Ньютон-Рафсон вы также сами программировали?
>> Да, сама
Это не самое лучшее решение. Я согласен, что полезно понять как что
происходит. Но дальше есть огромное количество мелких деталей.
Соответственно, с большой долей вероятности код из библиотеки будет
работать лучше. Зачем изобретать велосипед?
Вот например некоторые подпрограммы, которые использовал я при
разработке своей библиотеке
http://evgenii.rudnyi.ru/soft/tdlib00+/lib/toms/
Для минимизации была например lbfgsb.f.
Здесь только важно сделать свой код таким образом, чтобы можно было
легко менять численные подпрограммы. Тогда можно быстро пробовать разные
и найти лучшую. Это будет эффективнее, чем разрабатывать метод
оптимизации самой.
Наташка Перевощикова
Для нахождения начальных значений для метода Ньютона-Рафсона, я хочу использовать метод Нелдера — Мида (метод деформируемого многогранника и симплекс-метод). Как вы думаете, будеть ли он эффективен для вектора с размерностью больше, чем 14?Я читаю форумы по этому поводу и об этом методе отзываются положительно, но разные отзывы о его возможности искать минимум из-за размерности.
Если все ок с этим методом и , если Вы его использовали, то не могли бы Вы мне подсказать кайие ограничения накладывать на области, где интервал от 0 до 1?
Алгоритм взяла от сюда:http://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%9D%D0%B5%D0%BB%D0%B4%D0%B5%D1%80%D0%B0_%E2%80%94_%D0%9C%D0%B8%D0%B4%D0%B0
Evgenii Rudnyi
> симплекс-метод). Как вы думаете, будеть ли он эффективен для вектора с
> размерностью больше, чем 14?Я читаю форумы по этому поводу и об этом методе
> отзываются положительно, но разные отзывы о его возможности искать минимум
> из-за размерности.
Я никогда не использовал этот метод, поэтому ничего не могу сказать.
Надо попробовать. Только я бы посоветовал взять уже готовую подпрограмму.
Мне лично идея с выпуклой огибающей нравится больше. Я пробовал
использовать метод линейного программирования для нахождения начального
приближения. Это работало, но как это будет работать для
многокомпонентных систем сложно сказать. Надо также пробовать.
> Если все ок с этим методом и , если Вы его использовали, то не могли бы Вы
> мне подсказать кайие ограничения накладывать на области, где интервал от 0
> до 1?
Непростой вопрос. Если взять DONLP2 например, то можно в лоб задать это
как ограничения на параметры. Проблема будет тогда, когда мольная доля
при равновесии будет иметь маленькое значение, или будет близка к единице.
Можно сделать преобразование, см. например уравнение 6 в
http://evgenii.rudnyi.ru/doc/papers1/96calphad_bacuy.pdf
Мне кажется, что это также можно обобщить на многокомпонентную систему.
В этом случае область значений будет уже от минус до плюс бесконечности.
Проблема однако, что когда параметр будет слишком велик по абсолютной
величине (что соответствует очень маленькой мольной доле или значению
очень близком к единице), то все равно возникают численные проблемы.
Функция отклика перестает практически зависеть от параметра, что не
помогает успешному нахождению минимума.
Еще одно решение, которое я могу предложить, это замена логарифма. См.
ур. 3.9 в
http://evgenii.rudnyi.ru/soft/tdlib00+/doc/00tdlib.pdf
В этом случае мы позволяем мольной доле принимать отрицательные
значения, или значения больше единице, то есть область значений также
будет формально от минус до плюс бесконечности. При этом энергия Гиббса
ведет себя достаточно разумно, то есть минимум остается там где и надо,
и численно похоже все также неплохо. Однако, это идея осталась на
начальном уровне. Надо также пробовать.
Наташка Перевощикова
использовать метод линейного программирования для нахождения начального
приближения. Это работало, но как это будет работать для
многокомпонентных систем сложно сказать. Надо также пробовать.
Наташка Перевощикова
Evgenii Rudnyi
> использовать метод линейного программирования для нахождения начального
> приближения. Это работало, но как это будет работать для
> многокомпонентных систем сложно сказать. Надо также пробовать.
>
>
> Не могли бы Вы конкретнее написать какой метод Вы использовали?
Вот мой код
http://evgenii.rudnyi.ru/soft/bacuy_eq/
В качестве подпрограммы линейного программирования я использовал IMSL.
Мне очень нравилось. Далее для линейного программирования я использовал
бесплатную подпрограмму Minit
http://evgenii.rudnyi.ru/soft/tdlib00+/lib/minit/
но должен сказать, что в IMSL было сделано добротнее.
Evgenii Rudnyi
> оболочкам. Говорит, что я должна сама написать ее. Знать бы еше алгоритм.
> Какой-то кошмар.
Здесь по идее надо сказать, что надо выбрать, либо вы разрабатываете
численные методы для выпуклых оболочек, либо решаете термодинамические
задачи.
Но лучше всего делать так. Возьмите например алгоритм из Qhull, там
должны быть статьи, покажите его. Посчитайте число строк в Qhull.
Соответственно будет не сложно оценить время разработки. Далее надо
составить план на соотвествующее число лет и предложить его.