- 01
- 02
- 03
- 04
- 05
- 06
- 07
- 08
- 09
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
#include<stdio.h>
#include<math.h>
int const n = 50, n1 = 40, m = n / n1;
double f1(double t, double x, double y, double z);
double f2(double t, double x, double y, double z);
double f3(double t, double x, double y, double z);
int main() {
int i, k, ll;
double k11, k12, k13, k21, k22, k23, k31, k32, k33, k14, k24, k34;
double x, y, z, max;
double res[12][n1 + 1], h, a, b, t, pi, xn, yn, zn;
max = 1000.0; x = 3.0;
y = -2.0; z = -3.0;
a = 0.0; b = 1.0;
h = (b - a) / double(n);
k = 0; t = a;
res[0][0] = t; res[1][0] = x;
res[2][0] = y;
res[3][0] = z;
for (i = 1; i <= n; i++) {
k11 = f1(t, x, y, z);
k21 = f2(t, x, y, z);
k31 = f3(t, x, y, z);
k12 = f1(t + h / 2.0, x + h*k11 / 2.0, y + h*k21 / 2.0, z + h*k31 / 2.0);
k22 = f2(t + h / 2.0, x + h*k11 / 2.0, y + h*k21 / 2.0, z + h*k31 / 2.0);
k32 = f3(t + h / 2.0, x + h*k11 / 2.0, y + h*k21 / 2.0, z + h*k31 / 2.0);
k13 = f1(t + h / 2.0, x + h*k12 / 2.0, y + h*k22 / 2.0, z + h*k32 / 2.0);
k23 = f2(t + h / 2.0, x + h*k12 / 2.0, y + h*k22 / 2.0, z + h*k32 / 2.0);
k33 = f3(t + h / 2.0, x + h*k12 / 2.0, y + h*k22 / 2.0, z + h*k32 / 2.0);
k14 = f1(t + h, x + h*k13, y + h*k23, z + h*k33);
k24 = f2(t + h, x + h*k13, y + h*k23, z + h*k33);
k34 = f3(t + h, x + h*k13, y + h*k23, z + h*k33);
x = x + h*(k11 + 2.0*(k12 + k13) + k14) / 6.0;
y = y + h*(k21 + 2.0*(k22 + k23) + k24) / 6.0;
z = z + h*(k31 + 2.0*(k32 + k33) + k34) / 6.0;
t = t + h;
if (i%m == 0) {
k = k + 1;
res[0][k] = t;
res[1][k] = x;
res[2][k] = y;
res[3][k] = z;
res[4][k] = exp(t) + exp(2.0*t) + exp(-t);
res[5][k] = exp(t) - 3.0*exp(-t);
res[6][k] = exp(t) + exp(2.0*t) - 5.0*exp(-t);
res[7][k] = res[4][k] - res[1][k];
res[8][k] = res[5][k] - res[2][k];
res[9][k] = res[6][k] - res[3][k];
} if (res[7][k] < 0.0) {
res[7][k] = -res[7][k];
}
else if (res[8][k] < 0.0) {
res[8][k] = -res[8][k];
}
else if (res[9][k] < 0.0) {
res[9][k] = -res[9][k];
}
res[10][k] = res[7][k] + res[8][k] + res[9][k];
}
ll = k;
printf("k=%d\n", ll);
for (i = 0; i <= ll; i++) {
if (res[10][i] < max) max = res[10][i];
} for (i = 0; i <= n1; i++) {
printf("t=%.16lf\n", res[0][i]);
printf("x=%.16lf x(exact)=%.16lf delta=%.16lf\n", res[1][k], res[4][k], res[1][k] - res[4][k]);
printf("y=%.16lf y(exact)=%.16lf delta=%.16lf \n", res[2][k], res[5][k], res[2][k] - res[5][k]);
printf("z=%.16lf z(exact)=%.16lf delta=%.16lf \n", res[3][k], res[6][k], res[3][k] - res[6][k]);
}
printf("result:\n");
printf("t=%.16lf\n", res[0][n1]);
printf("x=%.16lf x(exact)=%.16lf delta=%.16lf\n", res[1][n1], res[4][n1], res[1][n1] - res[4][n1]);
printf("y=%.16lf y(exact)=%.16lf delta=%.16lf \n", res[2][n1], res[5][n1], res[2][n1] - res[5][n1]);
printf("z=%.16lf z(exact)=%.16lf delta=%.16lf \n", res[3][n1], res[6][n1], res[3][n1] - res[6][n1]);
printf("max norma|x|1=%.16lf \n", max / 3.0);
system("pause");
return 0;
}
double f1(double t, double x, double y, double z)
{
return x + z - y;
}
double f2(double t, double x, double y, double z)
{
return x + y - z;
}
double f3(double t, double x, double y, double z)
{
return 2.0*x - y;
}
Кандидат физико-математических наук сделал методичку по предмету "Моделирование систем". Это он так описал алгоритм решения системы ОДУ методом Рунге-Кутты.
Follow us!