Новичок
Профиль
Группа: Участник
Сообщений: 24
Регистрация: 4.5.2009
Репутация: нет
Всего: нет
|
Здравствуйте, форумчане!  Пишу программу для решения систем дифференциальных уравнений методом ФДН (формула дифференцирования назад). Проблема в том, что уже который день мучаюсь, но процесс постоянно расходящийся. Объясню проблему на примере. Необходимо решить задачу Коши для такой системы ДУ : y_ = y + 5*z z_ = - y - 3*z где y_, z_ - производные. Отрезок интегрирования [0,4], начальные условия y(0)=1, z(0)=1. Обычные одношаговые или многошаговые методы решают без проблем (или у меня просто не возникло там проблем). Сейчас программа как-то неправильно работает, т.к. решение сразу начинает "уходить" от настоящего. | Код |
#include <iostream>
using namespace std;
#include <vector>
#include <cmath>
#include "WSystemNewton.h"
#include "GaussMainColumn.h"
vector<double> X;
vector<double> Y;
vector<double> Z;
double f1(double x, double y, double z, double y_, double z_)
{
return -y_ + y + 5 * z;
}
double f2(double x, double y, double z, double y_, double z_)
{
return -z_ - y - 3 * z;
}
//--------------------------------------------------------------------------------//
double calcYn(int p, char func, double h, int k)
{
double* A = new double[(p + 1) * (p + 1)];
double* B = new double[p + 1];
double* gamma = new double[p + 1];
double xk1 = X[k] + h;
// заполняем первую строку
for (int j = 0; j < p + 1; j++)
A[j] = 1;
// заполняем 2..p строки матрицы
for (int i = 1; i < p + 1; i++)
{
for (int j = 0; j < p + 1; j++)
{
A[i * (p + 1) + j] = pow((xk1 - X[k - j]) / (h), i);
}
}
// заполняем вектор свободных членов
B[0] = 1;
for (int i = 1; i < p + 1; i++)
B[i] = 0;
GaussMainColumn gmc(p + 1, A, B);
// gmc.printMatrix();
gmc.setDebugMode(false);
gmc.process();
// gmc.printMatrix();
gamma = gmc.getX();
// for(int i = 0; i < p + 1; i++)
// cout << "gamma[i] = " << gamma[i] << endl;
// cout << "-------------------" << endl;
double yk1n = 0;
double gammaI = 0;
double yk1_i = 0;
if (func == 'y')
for (int i = 1; i <= p + 1; i++)
{
gammaI = gamma[i - 1];
yk1_i = Y[k + 1 - i];
yk1n += gammaI * yk1_i;
}
if (func == 'z')
for (int i = 1; i <= p + 1; i++)
{
gammaI = gamma[i - 1];
yk1_i = Z[k + 1 - i];
yk1n += gammaI * yk1_i;
}
// cout << "yk1n = " << yk1n << endl;
delete[] A;
delete[] B;
return yk1n;
}
//--------------------------------------------------------------------------------//
double calcY_(int p, char func, double h, int k, double yk1)
{
double* A = new double[(p + 1) * (p + 1)];
double* B = new double[p + 1];
double* alpha = new double[p + 1];
double xk1 = X[k] + h;
// заполняем первую строку
for (int j = 0; j < p + 1; j++)
A[j] = 1;
// заполняем 2..p строки матрицы
for (int i = 1; i < p + 1; i++)
{
for (int j = 0; j < p + 1; j++)
{
if (j != 0)
A[i * (p + 1) + j] = pow((xk1 - X[k - j + 1]) / (h), i);
else
A[i * (p + 1) + j] = 0;
}
}
// заполняем вектор свободных членов
for (int i = 1; i < p + 1; i++)
B[i] = 0;
B[1] = 1;
GaussMainColumn gmc(p + 1, A, B);
// gmc.printMatrix();
gmc.setDebugMode(false);
gmc.process();
// gmc.printMatrix();
for (int i = 0; i < p + 1; i++)
alpha[i] = gmc.getX()[i];
// for(int i = 0; i < p + 1; i++)
// cout << "alpha[i] = " << alpha[i] << endl;
// cout << "-------------------" << endl;
double yk1_ = 0;
double Yk1_i = 0;
if (func == 'y')
for (int i = 0; i <= p; i++)
if (i != 0)
{
Yk1_i = Y[k + 1 - i];
yk1_ += alpha[i] * Yk1_i;
}
else
yk1_ += alpha[i] * yk1;
if (func == 'z')
for (int i = 0; i <= p; i++)
if (i != 0)
yk1_ += alpha[i] * Z[k + 1 - i];
else
yk1_ += alpha[i] * yk1;
yk1_ *= -1.0 / h;
// out << "yk1_ = " << yk1_ << endl;
delete[] A;
delete[] B;
delete[] alpha;
return yk1_;
}
int main()
{
//--------------------------------------------------//
// Инициализация
double a = 0;
double b = 4;
double x0 = a;
double y0 = 1;
double z0 = 1;
double h = (b - a) / 1000;
// минимальный и максимальный шаги интегрирования
double hmin = (b - a) / 1000000.0;
double hmax = (b - a) / 4.0;
// минимальный и максимальный порядки методов
double pmin = 2;
double pmax = 6;
double p = 3;
WSystemNewton sys;
sys.f1 = f1;
sys.f2 = f2;
double xk = x0, yk = y0, zk = z0;
double xk1 = xk + h;
double zk1 = 0, yk1 = 0;
// коефициенты k для метода Рунге-Кутта 3 порядка
double k1y = 0, k1z = 0, k2y = 0, k2z = 0, k3y = 0, k3z = 0;
X.push_back(xk);
Y.push_back(yk);
Z.push_back(zk);
int k = 0;
cout << "xk1 = " << xk << "\tyk1 = " << yk << "\tzk1 = " << zk << endl;
//--------------------------------------------------//
// определение начальных значений для стартования
for (int i = 0; i < p; i++)
{
xk1 = xk + h;
// считаем коефициент k1 для всех уравнений
k1y = h * f1(xk, yk, zk, 0, 0);
k1z = h * f2(xk, yk, zk, 0, 0);
// считаем коефициент k2 для всех уравнений
k2y = h * f1(xk + h / 3.0, yk + k1y / 3.0, zk + k1z / 3.0, 0, 0);
k2z = h * f2(xk + h / 3.0, yk + k1y / 3.0, zk + k1z / 3.0, 0, 0);
// считаем коефициент k3 для всех уравнений
k3y = h * f1(xk + (2 / 3.0) * h, yk + (2 / 3.0) * k2y, zk + (2 / 3.0)
* k2z, 0, 0);
k3z = h * f2(xk + (2 / 3.0) * h, yk + (2 / 3.0) * k2y, zk + (2 / 3.0)
* k2z, 0, 0);
// искомые функции
yk1 = yk + (1 / 4.0) * (k1y + 3 * k3y);
zk1 = zk + (1 / 4.0) * (k1z + 3 * k3z);
xk += h;
yk = yk1;
zk = zk1;
// запоминаем значения решения
X.push_back(xk1);
Y.push_back(yk1);
Z.push_back(zk1);
cout << "xk1 = " << X[k + 1] << "\tyk1 = " << Y[k + 1] << "\tzk1 = "
<< Z[k + 1] << endl;
k++;
}
double yk1n = 0;
double zk1n = 0;
double yk1_ = 0;
double zk1_ = 0;
double alpha = 0;
xk = xk1;
while (xk1 + h <= b + 0.01)
{
xk1 = xk + h;
yk1n = calcYn(p, 'y', h, k);
zk1n = calcYn(p, 'z', h, k);
yk1_ = calcY_(p, 'y', h, k, yk1n);
zk1_ = calcY_(p, 'z', h, k, zk1n);
cout << "yk1n = " << yk1n << "\tzk1n = " << zk1n << endl;
cout << "yk1_ = " << yk1_ << "\tzk1_ = " << zk1_ << endl;
sys.setFirstCondition(xk1, yk1n, zk1n, yk1_, zk1_);
sys.calculate();
yk1 = sys.getY();
zk1 = sys.getZ();
xk = xk1;
X.push_back(xk1);
Y.push_back(yk1);
Z.push_back(zk1);
cout << "xk1 = " << X[k + 1] << "\tyk1 = " << Y[k + 1] << "\tzk1 = "
<< Z[k + 1] << endl;
cout << "------------------------------" << k << "------------------------------------" << endl;
k++;
}
return 0;
}
|
Если нужна дополнительная информация - спрашивайте. Очень нужна помощь. Спасибо! Это сообщение отредактировал(а) w1zard350 - 16.4.2010, 21:53 |
Формула дифференцирования назад, проблема сходимости 
Загрузка ...
