Численное решение задачи Коши, результат выводится неправильно - C (СИ)
Формулировка задачи:
не понимаю в чем проблема
задание было такое
Построить алгоритмы численного решения задачи Коши явным, неявным методом
Эйлера и методом Рунге-Кутта с автоматическим выбором шага интегрирования для
достижения заданной локальной относительной погрешности либо локальной
абсолютной погрешности . Апостериорную оценку погрешности для явного метода
Эйлера провести по правилу Рунге либо по оценке производной, входящей в
выражение для погрешности.
2. Составить рабочую программу с использованием универсальных функций численного
решения задачи Коши. Минимальный набор параметров функции решения задачи
Коши должен включать: имя функции, вычисляющей производную неизвестной
функции (правую часть дифференциального уравнения); интервал времени, на
котором нужно найти решение; начальное условие; погрешность решения; массив
временных отсчетов, в которых найдено решение; массив решений; максимальное
число временных отсчетов; количество фактически полученных временных отсчетов,
в которых найдены решения с заданной погрешностью.
#define _USE_MATH_DEFINES #include <stdio.h> #include <stdlib.h> #include <math.h> int a1; float TEC(float t, float U) { float z, e, E = 5, id, i0 = 3e-9, f = 1000, R = 5000, c = 1e-6, m = 1.7, fi = 0.026; a1++; e = E*sin(2 * (M_PI)*f*t + (M_PI) / 4); id = i0*(exp((e - U) / (fi*m)) - 1); z = id / c - U / (R*c); return z; } int RungeKutt(float f(float x, float y), float *t, float *y, float A, float B, float ny, float E, int Nmax) { int i, j, n = 0; float q = 2, R, y1, y2, y3, t1, t2, t3, h, a[6] = { 0, (1. / 4), (3. / 8), (12. / 13), 1., (1. / 2) }, c[6] = { 16. / 135, 0, 6656. / 12825, 28561. / 56430, -9. / 50, 2. / 55 }, c1[6] = { 25. / 216, 0, 1408. / 2565, 2197. / 4104, -1. / 5, 0 }, b[5][5] = { { 1. / 4, 0, 0, 0, 0 }, { 3. / 32, 9. / 32, 0, 0, 0 }, { 1932. / 2197, -7200. / 2197, 7296. / 2197, 0, 0 }, { 439. / 216, -8., 3680. / 513, -845. / 4104, 0 }, { -8. / 27, 2., -3544. / 2565., 1859. / 4104, -11. / 40 } }, k[6]; h = B - A; t1 = A; y1 = ny; do { y3 = y1; t3 = t1; for (i = 0; i<6; i++) { t2 = t1; y2 = y1; t2 += a[i] * h; for (j = 0; j<i; j++) { y2 += h*b[i - 1][j] * k[j]; } k[i] = f(t2, y2); } for (i = 0; i<6; i++) y1 += h*c[i] * k[i]; t1 += h; R = 0; for (i = 0; i<6; i++) R += (c[i] - c1[i])*k[i]; R = h*R; if (fabs(R) <= E*fabs(y1)) { y[n] = y1; t[n] = t1; if (fabs(R)<(E*fabs(y1) / 64)) h = h*q; if ((t1 + h)>B) h = B - t1; n++; } else { h = h / q; y1 = y3; t1 = t3; } } while ((t1<B) && (n<Nmax)); return n; } int yavniyEuler(float f(float x, float z), float *t, float *y, float a, float b, float ny, float E, float D, int Nmax) { int i, n = 0, q = 2; float h, y1, t1, y2, y3, ny1, nt1, d, k = 0, R; h = b - a; t1 = a; y1 = ny; do { ny1 = y1; nt1 = t1; if (k == 0) { d = f(t1, y1); y2 = y1 + h*d; } k = 1; for (i = 0; i<q; i++) { if (i == 0) { y1 = y1 + (h / q)*d; y3 = y1; } else y1 = y1 + (h / q)*f(t1, y1); t1 += (h / q); } R = E*fabs(y1); if (y2 == 0) { R = R + D; } if (fabs(y1 - y2) <= (R)) { k = 0; y[n] = y1; t[n] = t1; if (fabs(y1 - y2)<E*fabs(y1) / 4) h = q*h; n++; } else { h = h / q; y2 = y3; y1 = ny1; t1 = nt1; } } while (t1<b && n<Nmax); return (n); } int neyavniyEuler(float f(float x, float z), float *t, float *y, float a, float b, float ny, float E, int Nmax) { int n = 0, q = 2; float h, h1, d, em = 0.000000119, y1, t1, y2, t2, y3, y4, f1, f2, v, R, R1; h = b - a; y1 = ny; t1 = a; y2 = a; d = f(a, ny); do { y4 = y1 + h*d; t2 = t1 + h; y3 = y2; if (y2 != 0) h1 = y2*sqrt(em); else h1 = sqrt(em); f1 = y2 - y1 - h*f(t2, y2); f2 = (y2 + h1) - y1 - h*f(t2, (y2 + h1)); y2 = y2 - f1*h1 / (f2 - f1); while (fabs(y2 - y3)>E*fabs(y2)) { f2 = f1; f1 = y2 - y1 - h*f(t2, y2); v = y2; y2 = y2 - (y2 - y3)*f1 / (f1 - f2); y3 = v; } R = (y4 - y2) / 2; if (fabs(R) <= E*fabs(y2)) { y1 = y2; t1 = t2; y[n] = y2; t[n] = t2; d = f(t1, y1); if (fabs(R)<E*fabs(y2) / 4) h = h*q; if ((t2 + h)>b) h = b - t2; n++; } else { h = h / q; } } while (t1<b && n<Nmax); return n; } int main(void) { int j, i, Nmax = 10000; float *t, *y; t = malloc(Nmax*sizeof(float)); y = malloc(Nmax*sizeof(float)); a1 = 0; i = yavniyEuler(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 0.01, 10000); printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1); a1 = 0; i = neyavniyEuler(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 0.01, 10000); printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1); a1 = 0; i = RungeKutt(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 10000); printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1); }
Решение задачи: «Численное решение задачи Коши, результат выводится неправильно»
textual
Листинг программы
float f(float x, float y)
ИИ поможет Вам:
- решить любую задачу по программированию
- объяснить код
- расставить комментарии в коде
- и т.д