Решение уравнения в точках 0,25; 0,5; 0,75; Найти точное решение оду

Задание по лабораторной работе:
Вычислить с заданной точностью (0.001) решение уравнения в точках 0,25; 0,5; 0,75; 1.

Найти точное решение ОДУ.

Сравнить приближенное решение с точным.
Задача Коши, согласно заданию лабораторной работы, имеет вид

Найдем точное решение задачи Коши. Поскольку мы имеем линейное неоднородное уравнение второго порядка, найдем сначала общее решение линейного однородного уравнения, решая характеристическое уравнение

. Его решениями являются 0 и 1, так что общее решение линейного однородного уравнения имеет вид

. Чтобы найти частное решение неоднородного уравнения, будем искать его в виде

, так как у характеристического уравнения два вещественных корня, один из которых 0, а правая часть – в виде полинома.

Получим, что частное решение неоднородного уравнения имеет вид, а общее решение неоднородного уравнения имеет вид. Подставляя начальные значения, находим точное решение поставленной задачи Коши:

Теперь приступим к численному решению задачи Коши.

Численное решение методом Рунге-Кутта в нашем случае облегчается тем, что фактически требуется задание только одной функции, поскольку

Для решения задачи воспользуемся языком программирования C++, а функцию для решения методом Рунге-Кутта сделаем такой, что она будет не просто вычислять новые значения y и z, но и записывать их в старые переменные, что позволит легко использовать ее в циклах. Для обобщенности решения будем передавать ей в качестве параметра указатель на указанную выше функцию

.
Вот как выглядит код этой функции для вычисления очередного шага методом Рунге-Кутта:

double RK(double& x, double& y, double& ys, double h,

double (*f)(double,double,double))
{
double k1,k2,k3,k4;
double m1,m2,m3,m4;
k1 = h*f(x,y,ys);
m1 = h*(ys);
k2 = h*f(x+h/2.,y+m1/2.,m2 = ys+k1/2.);
m2*=h;
k3 = h*f(x+h/2.,y+m2/2., m3 = ys+k2/2.);
m3*=h;
k4 = h*f(x+=h, y+m3, m4 = ys+k3);
m4*=h;
ys+=(k1+2.*(k2+k3)+k4)/6.;
return y+=(m1+2.*(m2+m3)+m4)/6.;
};
Полностью автоматизируем решение нашей задачи, для чего создадим два массива — для решений в интересующих нас точках с шагом h и с шагом h/2, и организуем вычисления в цикле, со все уменьшающимся шагом. После каждого вычисления будем проверять, достигнута ли требуемая точность, и если нет — то цикл будет повторен с уменьшенным значением (перед этим следует переписать значения из массива для половинного шага в первый массив).

Таким образом, расчетная часть программы имеет вид:

double Y(double x, double y, double ys)
{
return 3+ys;
}

int main(int argc, const char * argv[])

{
double eps = 0.0001;
double yh[4], zh[4], yh2[4], zh2[4];
double h = 0.25; // Начальное значение шага
int N = 4; // Количество шагов
int n = 1; // Записывать каждое n-е значение

// Заполняем таблицу для шага h

double x = 0, y = 6, z = 2;
for(int i = 1; i <= N; ++i)
{
RK(x,y,z,h,Y);
if (i%n == 0) {
yh[i/n-1] = y;
zh[i/n-1] = z;
}
}

// Цикл, пока не достигнем точности.

for(;;)
{
// Уменьшаем шаг
N *= 2; n *= 2; h /= 2.0;
x = 0; y = 6; z = 2;
for(int i = 1; i <= N; ++i)
{
RK(x,y,z,h,Y);
if (i%n == 0) {
yh2[i/n-1] = y;
zh2[i/n-1] = z;
}
}
// Проверяем точность
bool ok = true;
for(int i = 0; i < 4; ++i)
{
if (fabs(yh[i] - yh2[i])/15 > eps)

{ ok = false; break; }

if (fabs(zh[i] - zh2[i])/15 > eps) { ok = false; break; }
}
if (ok == false)
{
// Точности не хватает!
// Переписываем массив и считаем с новым шагом
for(int i = 0; i < 4; ++i)
{
yh[i] = yh2[i];
zh[i] = zh2[i];
}
continue;
}
break;
}

// Выводим таблицу значений для данного шага

x = 0; y = 6; z = 2;
printf(" x y y' y точное"
" Погрешность y' точное Погрешность\n");
char * const fmt = "%.6lf %.6lf %.6lf "
"%.6lf %10.3lg %.6lf %10.3lg\n";
printf(fmt,x,y,z,
5*exp(x)+1-3*x,5*exp(x)+1-3*x-y,

5*exp(x)-3,5*exp(x)-3-z);

for(int i = 0; i < N; ++i)
{
RK(x,y,z,h,Y);
printf(fmt,x,y,z,
5*exp(x)+1-3*x,5*exp(x)+1-3*x-y,

5*exp(x)-3,5*exp(x)-3-z);

}

}
Результат работы программы:

x y y' y точное Погрешность y' точное Погрешность

0.000000 6.000000 2.000000 6.000000 0 2.000000 0

0.125000 6.290741 2.665741 6.290742 1.3e-006 2.665742 1.3e-006

0.250000 6.670124 3.420124 6.670127 2.94e-006 3.420127 2.94e-006

0.375000 7.149952 4.274952 7.149957 5e-006 4.274957 5e-006

0.500000 7.743599 5.243599 7.743606 7.56e-006 5.243606 7.56e-006

0.625000 8.466219 6.341219 8.466230 1.07e-005 6.341230 1.07e-005

0.750000 9.334986 7.584986 9.335000 1.46e-005 7.585000 1.46e-005

0.875000 10.369357 8.994357 10.369376 1.92e-005 8.994376 1.92e-005

1.000000 11.591384 10.591384 11.591409 2.49e-005 10.591409 2.49e-005

Как видно из расчета, при шаге 0.125 погрешность расчета существенно меньше указанной в задании 0.0001.