где
и - параметры спектральной плотности,
, , и - коэффициенты уравнения,
и начальными условиями:
и временем моделирования 120 сек, причем относительная погрешность среднеквадратического отклонения ,
если:
а) случайное воздействие имеет спектральную плотность ;
б) если случайное воздействие X(t) является белым шумом.
Моделирование выполняется с целью вычисления количества ординат случайного процесса y(t), которые выходят за уровень
Выполним математическое моделирование непрерывно-стохастической системы.
Будем использовать нелинейное стохастическое уравнение 2-го порядка , (1)
где - случайный процесс.
Для реализации математической модели в случаях:
а) случайное воздействие имеет спектральную плотность , (2)
- круговая частота;
- коэффициент затухания корреляционной функции;
- средняя частота корреляционной функции.
а) если случайный процесс имеет спектральную плотность.
Белый шум - стационарный случайный процесс с нулевым математическим ожиданием и корреляционной функцией, равной дельта-функции.
Моделирование белого шума осуществляется по следующей формуле:
, (3)
-независимая случайная величина с нормальным законом распределения с mx=0 и Dx=1,
No - коэффициент интенсивности белого шума или высота спектральной плотности.
Моделирование случайного воздействия со спектральной плотностью осуществляется стохастическим дифференциальным уравнением второго порядка
; (4)
в систему уравнений 1-ого порядка, для этого введем специальные переменные:
(5)
В результате получим следующую систему 1-го порядка:
(6)
Применяем к каждому уравнению метод Эйлера
(7)
получим следующую численную модель:
(8)
В случае а) когда случайное воздействие – белый шум, аналогично, математическая модель будет иметь вид:
(9)
При моделировании непрерывной стохастической модели следует выполнить такие действия:
1) Подбор коэффициента интенсивности белого шума (его мы осуществим с помощью табуляции функции
,
ее максимальное значение и будет требуемым шагом);
2) разработать датчик случайных чисел с нормальным законом распределения.
Для этого необходимо:
- сгенерировать два случайных числа с равномерным законом распределения, 1-ое число , а второе число
(Рисунок 1);
- сравнить, если V1>f(V1), то все числа отбрасываются и генерация повторяется заново, иначе меньшее число принимается как верное;
3) выбрать произвольный шаг табулирования;
4) получить значения по системам уравнений (8),(9);
5) проверить сходимость - проверка выполняется среднеквадратично по формуле
, (10)
Если погрешность среднеквадратичного отклонения менее или равна 0.05, то полученные значения считаются решением, иначе необходимо уменьшить шаг в 2 раза и повторить итерацию.
Причем в случае, где X(t)- белый шум обеспечиваем сходимость только по x1 (8); а в случае, где случайное воздействие имеет спектральную плотность (2), сходимость обеспечиваем и по x1 и по x3.
На основе выбранной численной модели была разработана программа по моделирования системы.
Алгоритм работы программы следующий:
- находится коэффициент интенсивности белого шума No, для этого функция табулируется , в диапазоне (1;120) с шагом 0,1
Первая часть задачи, где m(t) белый шум:
- применяется генератор случайных чисел с нормальным распределением;
- выбирается произвольный шаг;
- получаются зависимости y(t) от t и y’(t);
- выполняется контроль среднеквадратического отклонения
по формуле
-если среднеквадратического отклонения менее, либо равно 0.05 то полученные зависимости считаются решением, иначе шаг табулирования уменьшается в два раза.
Решение второй части задачи, где х(t) заданная функция, выполняется по выше описанному алгоритму лишь с той разницей, что контроль среднеквадратического отклонения ведется не только по x1, но и по x3. (из формулы (6 ) ). Полученный результат выводится в текстовый файл.
После завершения работы программы были получены необходимые точечные оценки дифференциального стохастического уравнения.
Результаты представлены ниже на рисунках 1-6.
Программа приведена в приложении А.
Случайный процесс является белым шумом:
Рисунок 1- Зависимость y от t
Рисунок 2 - Зависимость y’ от t
Случайное воздействие на систему- заданная функция:
Рисунок 3 – Зависимость y от t
Рисунок 4 – Зависимость z от t
Была выполнена работа по моделированию состояния системы непрерывно-стохастической модели на ЭВМ, состояние которой описывается стохастическим дифференциальным уравнением ,
со следующими параметрами:
, , и -коэффициенты уравнения,
и временем моделирования 120 сек, относительная погрешность среднеквадратического отклонения,
б) если случайный процесс является белым шумом.
В данной работе:
Ø выбрали метод моделирования стохастической дифференциальной системы;
Ø построили численную модель состояния системы;
Ø выполнили моделирование по построенной численной модели;
Ø оценили выброс случайной величины за уровень ;
Ø Выполнили проверку датчика сл.чис.с помощью критерия Хи квадрат.
(Приложение Б)
Моделирование выполнялось с целью вычисления количества ординат случайного процесса y(t), которые выходят за уровень и подсчет количества выхода значений за этот уровень – ни одно значение не вышло за уровень.
1. Томашевский В. М., Жданова В. Г., Жолдаков О.О.. Вирішення практичних завдань методами комп’ютерного моделювання: Навч. посібник.- К.:”Корнійчук”,2001.-268с.
2. Статистические методы для ЭВМ/ Под ред. К.Энслейна: Пер. с англ. / Под ред. М.Б.Малютова.- М.:Наука. Гл.ред. физ. Мат., лит. 1986.-464с.
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int main(){
int const k=1000;
double t,y,z;
int i,j;
int n=120;
int n0=1;
double w=1;
double b1=0.5;
double b2=1.5;
double c1=1.2;
double c3=-1.5;
double M=0.03;
double h=0.1;
t=0;
double t1=0;
z=0;
double z1=0;
y=0.15;
double y1=0.15;
double y_max,y_rez,eps,eps1;
double mas1[1200];
double mas[1200];
double e;
FILE *stream;
printf(" t | y | z \n");
/// fprintf(stream, " t | y | z \n");
i=0;
/* open a file for update */
stream = fopen("DUMMY.FIL", "w+");
while (t<120)
{
double j1, j2, r1, r2, s;
j1 = -1.0 + 2.0*rand()/((double)RAND_MAX - 1.0);
j2 = -1.0 + 2.0*rand()/((double)RAND_MAX - 1.0);
s = j1*j1 + j2*j2;
if (s < 1)
r1 = j1 * sqrt(-2*log(s)/s);
r2 = j2 * sqrt(-2*log(s)/s);
};
e=r1;
t=t+h;
y=y+h*z;
z=z+h*(e*pow((n0/h),1/2)-b1*z-b2*fabs(z)*z-c1*y-c3*pow(y,3));
mas[i]=y;
printf("%f | %f | %f \n",t,y,z);
/* write some data to the file */
fprintf(stream, "%f | %f | %f \n",t,y,z);
t1=t1+h/2;
y1=y+((h/2)*z1);
z1=z1+((h/2)*(e*pow((n0/h),1/2)-b1*z-b2*fabs(z)*z-c1*y-c3*pow(y,3)));
mas1[i]=y1;
i=i+1;
}
/* close the file */
// fclose(stream);
y_rez=0;
y_max=mas1[0];
for (i=0; i<=119;i++)
y_rez=y_rez+(pow((mas1[i]-mas[i]),2));
if (mas1[i]>y_max)
y_max=mas1[i];
eps=pow(y_rez/n,0.5)/fabs(y_max);
printf("%f epsilon for our operations\n",eps);
fprintf(stream, "%f epsilon for our operations\n",eps);
fclose(stream);
return 0;
Страницы: 1, 2
