Гуд бай, Blogger.com!

Этот так называемый веблог переехал на

http://structural.wordpress.com/

со всеми сообщениями. Ура!

Matlab: статистические штучки для получения реализаций случайной величины

Это первое сообщение, отосланное через ScribeFire, аддон к браузеру FireFox.
Так-с. Исходные данные. Предполагается, что читатель знаком с:

• методом Монте-Карло (вкратце: случайные параметры -> много-много расчетов посредством произвольной расчетной модели -> вероятностный результат, надежность и прочие радости)
• и с тем, что такое функция распределения и функция плотности случайной величины.

Предположим, что у нас есть случайная прочность R с параметрами m=350, s=30. Это выглядит так:

Функции распределения и функции плотности в Матлабе соответствуют команды

cdf('norm', x, Mean1, StdDev1)
pdf('norm', x, Mean1, StdDev1)

Команда cdf(...) отвечает на вопрос "какова вероятность того, что реализация случайной переменной будет меньше, чем x?"
Но для задачи вероятностного расчета методом Монте-Карло нам нужно средство генерации множества реализаций нашей случайной прочности.
Для равномерного распределения такое средство одинаково в любом языке программирования и называется RND. И в Екселе тоже есть.
Для любого другого распределения для генерации случайных чисел нам нужна функция, обратная CDF, отвечающая на вопрос "какое число соответствует вероятности непревышения такой-то?". Такой способ получения случайных чисел называется методом обратного преобразования.
Откуда нам взять такого зверя? На помощь приходит замечательный статистический пакет Stixbox, о котором я обещал написать :

>> MonteCarloData1 = rnorm(1000, 350, 30);
>> Mean1 = mean(MonteCarloData1)
Mean1 =
349.9133
>> StdDev1 = std(MonteCarloData1)
StdDev1 =
30.7880

Функция rnorm(количество_реализаций, среднее, стд_отклонение) выдает столько реализаций случайной величины с нормальным распределением, сколько нам заблагорассудится. С ее помощью и построена наложенная на рисунок гистограмма реализаций случайной величины (красненьким).
Аналогичные функции есть в этом пакете и для других известных распределений. Для задач вероятностного анализа расчетных моделей конструкций более чем достаточно. Спасибо за этот замечательный инструмент Андерсу Хольтсбергу.

Про 0,99865

Как показало былинное обсуждение «коэффициента спокойного сна» на dwg.ru, насчет критериев обеспеченности расчетных значений существует значительная неразбериха. Не говоря уже про вероятностный анализ.
Итак. Откуда берутся упомянутые в первой ссылке величины обеспеченности 0,99865 и 0,997. (Предполагается, что читатель знает, что такое среднее значение случайной величины и её стандартное отклонение.)

• 0,997 - вероятность того, что реализация случайной величины с нормальным распределением окажется в промежутке между значением "среднее минус три стандартных отклонения" и "среднее плюс три стандартных отклонения". Этот промежуток будет называться доверительным интервалом.
• 0,99865 - вероятность того, что эта реализация будет больше, чем значение "среднее минус три стандартных отклонения".

Первое можно получить с помощью так называемой функции ошибки (используем Matlab/Octave, как обычно):

>> erf(3/sqrt(2))
ans = 0.9973

Второе, как уже обсуждалось в более ранней записи,

>> 1 - cdf('norm', 1000-3*50, 1000, 50)
ans =   0.998650101968370

Как видим, 0,997 пришло к нам из статистики и еретического учения маркетинга сотоварищи. 0,99865, наоборот, есть правильная инженерная цифирь и, попросту говоря, является обеспеченностью (вероятностью), например, расчетной прочности бетона.
В нормах, понятно, все это отражено крайне мутными словами, и расчетное значение вычисляется на основе нормативного (с обеспеченностью 0,95) при помощи частных коэффициентов. Все это приводит к тому, что доценты нашего института не знают, что такое "обеспеченность" и в чем разница между коэффициентами перегрузки и коэффициентами безопасности (правильнее - частные коэффициенты надежности).

На картинке показана:
темно-красным - плотность распределения некоторой прочности с матожиданием 350, стандартным отклонением 30;
красным - гистограмма дискретных реализаций этой величины (1000 штук);
синим - расчетная и нормативная прочности с обеспеченностью 0,99865 и 0,95.

Про численные методы.

Про численные методы есть замечательные материалы и примеры для нескольких программных пакетов на сайте
http://numericalmethods.eng.usf.edu/

Matlab: количество отображаемых знаков

Этот пример считает обеспеченность величины, находящейся на расстоянии трех стандартных отклонений от среднего значения нормального распределения:

>> 1 - cdf('norm', 1000-3*50, 1000, 50)
ans =
    0.9987

Короче говоря, с инженерной точки зрения это обеспеченность расчетного значения случайной величины со средним значением 1000 и стандартным отклонением 50. Как видим, Матлаб съел все знаки после четвертого, хотя на самом деле посчитал точное число. Чтобы не путаться, можно отформатировать вывод Матлаба, чтобы он показывал число с точностью 'long':

>> format('long')
>> 1 - cdf('norm', 1000-3*50, 1000, 50)
ans =
   0.998650101968370

Вот она, наша знаменитая обеспеченность 0,99865.
В следующий раз будет
* про то, чем отличается 0,997 от 0,99865,
* про замечательный статистический аддон Stixbox,
* про инвертирование с его помощью функции распределения для генерации случайных чисел в методе Монте-Карло.

Короче, полезная в расчетах надежности статистика.
Но самое интересное впереди. Самое интересное будет про нечёткие множества.

Блог жив

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

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

Если вы это читаете, то вы наверняка с dwg.ru. Хе-хе.

Вывод предупреждений в ANSYS

Заглушить вывод предупреждений можно командой

/uis,msgpop,3

Отглушить обратно -

/uis,msgpop,0

Проблема управления памятью в Matlab

В Matlab есть одна неприятная штука -- Out of memory error. Фактор управления памятью приходится учитывать при операциях с большими матрицами, которые требуют немпрерывных массивов памяти.
Существует ряд способов решения проблем, связанных с выделением памяти.
Читать про них:
Memory Management Guide
Так, сегодня проблему с очень большой матрицей вектора перемещений в методе конечных элементов UX(i,j) размером 836x2007 мне удалось решить, используя очень простую запись в виде массиыва ячеек UX{i}(j) (массив ячеек хранит указатели на массивы, а не всю матрицу сразу. Просто? А вот я полдня с со своей программой мучился).
Следует отметить, что Matlab вываливался с ошибкой от матрицы, считываемой с диска с явно заданными элементами:

UX(1,1)=5.8;
UX(1,2)=5.4;
UX(1,2)=5.2;
...

Если вы создадите случайную матрицу rand(836,2007) прямо в командной строке, никакой ошибки не случится.

Про массивы ячеек (cell arrays) можно почитать очень дельную статью:
Cell Arrays and their contents

Как в Matlab соединить два текстовых файла

Пусть мы записываем что-нибудь в файл file1.txt и нам надо добавить туда содержимое файла file2.txt. Это выглядит следующим образом:

fid = fopen('file1.txt', 'wt'); %открыли для записи
fid2 = fopen('file2.txt', 'r'); %открыли для чтения
while 1
    tline = fgetl(fid2);
    disp(tline)
    fprintf(fid, '%s\n', tline);
    if feof(fid2), break, end
end
fclose(fid2);

Две каноничные страницы советов по ANSYS/APDL

Совершенно банальные для программиста ANSYS ссылки:

ansys.net Tips.

Peter Budgell's ANSYS Tips & Tricks page из лохматого 1999 года.

И Post Scriptum о графическом интерфейсе

Небольшая и полезная статья "Customizing GUI Behavior [PDF]" на английском языке с сайта ansys.net.

Доступ к средствам отрисовки в ANSYS

Через язык программирования (или правильнее называть язык команд?) APDL можно делать многие интересные вещи, которые нам пригодятся для осуществления хотя бы минимальных действий по визуализации результата.
Итак, цитата с форума dwg.ru:

Цитата:
Сообщение от Александр Бауск
я хочу, например, обозначить часть элементов (пластины, с известными номерами) одним цветом, часть -- другим. .

Ответ дорогого товарища The_Mercy_Seat:
Сначала селектируете элементы по нужному признаку, а потом назначаете им цвет U_M: PlotCtrls => Style => Colors => Picked Entity Colors. Потом выделяете все.

Итак, мораль: командой /COLOR можно по желанию расцвечивать наши элементы или другие объекты базы данных ANSYS.

XML_IO для MATLAB

Замечательная вещь для экспорта-импорта XML в среде MATLAB.
Написана полностью на языке MATLAB, очень удобно. Потом, если получится, отпишу о ее использовании в своем пакете интервального анализа.

Сама утилита
Небольшая демонстрация и руководство

Спасибо автору:
Jaroslaw Tuszynski
jaroslaw.w.tuszynski@saic.com

Кладезь для Ansys

Ну, вообще-то весь сайт www.ansys.net -- просто кладезь мудрости, даже по сравнению с роскошной пользовательской документацией ANSYS.
Недокументированные функции, справочник команд, советы, статьи по применению, в частности, очень меня беспокоящего SOLID65.

Ссылки на электронные библиотеки

"Библиотека Ихтика"

Книги на сайте зайберколбасы

Excel XML Toolbox

Продолжаем конструировать набор скриптов для обработки модели по схеме AutoCAD --> ANSYS --> Excel. Итак, Excel XML Toolbox для импорта данных из результатов ANSYS:

Download EXcel XML Toolbox

ANSYS. Вывод различных данных в файл XML

Недокументированная команда /XML и XMLOPT:
/XML,filename,ext,directory
This ANSYS command allows to write current load step results in XML format for use in Workbench Simulation.
Similar to other ANSYS commands, the 2nd through 4th arguments allow specification of the filename to save to. "Filename" defaults to the jobname, "Ext" defaults to "xml", and "Directory" defaults to the current working directory. Use of this command will write out geometry and results, as specified by the XMLOPT undocumented command.


Ссылки по теме:
http://ansys.net/collection/567