Метод сглаженных частиц (SPH). Задача Римана.

[zeher]

Здравствуйте. Проблема такая. Я написала программу для решения задачи Римана методом SPH. Но она не работает. И сложность в том что определить в чем ошибка (в формуле или в коде) я не могу. Прошу помощи. Если у кого то есть рабочая программа, поделитесь пжл, чтобы имела возможность сравнить промежуточные результаты и найти ошибку.
Система уравнений:
$ \rho_{i} = \sum\limits_{j}m_{j}\left(\vec{u}_{i}-\vec{u}_{j}\right)\vec{\nabla}_{i}W_{ij},$

$ \frac{dE_{i}}{dt}=\frac{1}{2}\sum\limits_{j}m_{j}\left(\frac{p_{j}}{\rho^{2}_{j}}+\frac{p_{j}}{\rho^{2}_{j}}\right)\left(\vec{u}_{i}-\vec{u}_{j}\right)\vec{\nabla}W_{ij},$


$ \frac{d\vec{u}_{i}}{dt}=-\sum\limits_{j}m_{j}\left(\frac{p_{j}}{\rho^{2}_{j}}+\frac{p_{j}}{\rho^{2}_{j}}\right)\vec{\nabla}W_{ij},$

$p=(1-\gamma)\rho E.$

Кратко описываю программу.
Имеется два класса. первый "частица", с параметрами - координата(одномерный случай), плотность, давление, энергия и скорость. Второй класс "метод sph", в нем все вычисления происходят.
Начальные значения параметров задаются в конструкторе. Остальные считываются с файла const.
Интегрирование проводится по координате, энергии и скорости. Методом эйлера.
Результат заносится в файл.
Самое интересное в том, что на первый шагах результат есть(не знаю правильный или нет), а потом значения неопределенные становятся.

Очень надеюсь на помощь. Я в полном отчаянии.
Вот ссылка на программу
http://rghost.ru/38307097