Решение параболических уравнений

n = ceil (l / h);

for (i = 0; i <= 119; i++)

{

for (j = 1; j <= 119; j++)

{

u[i] [j] = 0;

a[i] [j] = 0;

b[i] [j] = 0;

}

}

for (i = 0; i <= n; i++)

{

u[i] [0] = phi (i * h);

}

for (j = 0; j <= floor (T /tau); j++)

{

u[0] [j] = mu_1 (tau * j);

u[n] [j] = mu_2 (tau * j);

}

for (j = 0; j <= floor (T / tau); j

++)

{

a[1] [j + 1] = 1 / (2 + s);

for (i = 2; i <= n – 1; i++)

{

a[i] [j + 1] = 1 / (2 + s – a [i – 1] [j + 1]);

}

b[1] [j + 1] = mu_1 ((j + 1) * tau) + s * u[1] [j] + pow (h, 2) * f (h, (j + 1) * tau);

for (i = 2; i <= n – 1; i++)

{

b[i] [j + 1] = a [i – 1] [j + 1] + s * u[i] [j] + pow (h, 2) * f (i * h, (j + 1) * tau);

}

u[n] [j + 1] = mu_2 ((j + 1) * tau);

for (k = 1; k <= n – 1; k++)

{

u [n – k] [j + 1] = a [n – k] [j + 1] * (b [n – k] [j + 1] + u [n – k + 1] [j + 1]);

}

}

myfile = fopen («output.txt», «w+»);

plotter = fopen («3dplot.txt», «w+»);

fprintf (myfile, «Таблица значений функции u=u (x, t) в области D={0<=X<=%g, 0<=T<=%g}:\n», l, T);

printf («Значения функции u (x, t) в области D={0<=X<=%g, 0<=T<=%g}:\n\n», l, T);

for (j = 0; j <= floor (T / tau); j++)

{

for (i = 0; i <= n; i++)

{

printf («%.2f», u[i] [j]);

fprintf (myfile, «u(%g) (%g)=%g;\n», i * h, j * tau, u[i] [j]);

if (i < n && j < floor (T / tau))

{

fprintf (plotter, «[[%g, %g, %g], [%g, %g, %g], [%g, %g, %g], [%g, %g, %g]]», i * h, j * tau, u[i] [j], (i + 1) * h, j * tau, u [i + 1] [j], i * h, (j + 1) * tau, u[i] [j + 1], (i + 1) * h, (j + 1) * tau, u [i + 1] [j + 1]);

if (i >= n – 1 && j >= floor (T / tau) – 1)

{

}

else

{

fprintf (plotter,»,»);

}

}

}

printf («\n»);

}

fclose(myfile);

fclose(plotter);

printf («\nОсь x расположена горизонтально; ось t расположена вертикально и направлена вниз»);

printf («Шаг по оси x равен % g; шаг по оси t равен % g.\n», h, tau);

printf («\nДля выхода нажмите ENTER…»);

while (getch()!= 13);

}

// – //

float f (float x, float t)

{

return x * t;

}

// – //

float mu_1 (float t)

{

return 2.1 + t;

}

// – //

float mu_2 (float t)

{

return 3.2 * (t + 1 / 2.71828);

}

// – //

float phi (float x)

{

return (1.1 * pow (x, 2) + 2.1) * exp(-x);

}

 Скачать реферат