Уравнения в частных производных (стр. 2 )

Теперь, учитывая условие в (8.12), из (8.18), (8.19) выводим формулу для вычисления значений функции на первом слое:

. (8.20)

С учетом (8.13), окончательно получим для приближенных значений искомой функции на первом слое формулы

. (8.21)

Учитывая граничные условия (8.13) из (8.17) выводим формулы для вычисления значений на слоях :

Мы получили явные формулы решения разностной задачи.

Решение задачи называется устойчивым, если малым изменениям исходных данных отвечают малые изменения решения.

Приведем без доказательства (доказательство можно найти в [9]) следующий факт: для выполнения условия устойчивости разностной схемы (8.17) необходимо и достаточно, чтобы выполнялось условие Куранта cτ < h.

Сформулируем алгоритм решения задачи о колебаниях струны (8.11) — (8.13) с помощью явной разностной схемы (8.17).

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

1. Построить сеточную область

(8.21)

выбирая шаги h и τ так, чтобы выполнялось условие устойчивости (условие Куранта) cτ < h.

2. Вычислить значения ui,k искомой функции для k = 0, 1:

(8.22)

(8.23)

3. Вычислить значения ui,k для :

(8.24)

Пример 8.1. Решить модельную задачу

и проверить алгоритм (8.21) — (8.24). Точным решением задачи является функция .

Решение в Excel.

Выберем шаг по переменной x: h = 0,1; шаг по переменной t: τ = 0,05. Так как c = 1, то выполняется условие Куранта cτ < h. Поэтому можем применять расчетные формулы (8.22) — (8.24).

1. В ячейки B1:L1 запишем значения переменной x: 0; 0,1; …; 1.

2. В ячейки A2:AL22 запишем значения переменной t: 0; 0,05; …; 1.

3. В ячейку B2 введем формулу =SIN(1+B1), выделим B2 и протянем маркером заполнения вправо до ячейки L2. В этой строке мы получим ui,0.

4. В ячейку B3 введем формулу =SIN(1)*EXP(2*$A3), выделим B3 и протянем маркером заполнения вниз до ячейки B22. В этом столбце мы получим u0,k.

5. В ячейку L3 введем формулу =SIN(2)*EXP(2*$A3), выделим L3 и протянем маркером заполнения вниз до ячейки L22. В этом столбце мы получим un,k.

6. В ячейку С3 введем формулу

=SIN(1+C1)+2*SIN(1+C1)*0,05+0,05^2*((D2-2*C2+B2)/0,01+5*SIN(1+C1)*EXP(0))/2,

выделим С3 и протянем маркером заполнения вправо до ячейки K3. В этой строке мы получим ui,1.

7. В ячейку С4 введем формулу

=2*C3-C2+0,05^2*((B3-2*C3+D3)/0,01+5*SIN(1+C$1)*EXP(2*$A3)),

выделим С4 и протянем маркером заполнения вправо до ячейки K4. Затем выделим диапазон С4: K4 и протянем маркером заполнения вниз до 22-й строки.

В диапазоне B2:L22 получим значения сеточной функции ui,k. Чтобы проверить правильность решения, вычислим в 24-й строке значения точного решения задачи . Для этого введем в ячейку B24 формулу =SIN(1+B1)*EXP(2), затем протянем B24 маркером заполнения до L24. Значения в этой строке должны совпадать с соответствующими значениями строки B22:L22.

Вычислим в 25-й строке относительные погрешности. Введем в ячейку B25 формулу =ABS((B24-B22)/B24), выделим B25 и протянем маркером заполнения до L25. Мы увидим, что относительные погрешности меньше 0,00025.

На рис.8.3 представлен график решения задачи, построенный в программе Mathcad по таблице B2:L22.

Рис. 8.3

Составим функции в Excel для реализации алгоритма (8.22) — (8.24). Для данного примера имеем:

.

Выполним команду «Сервис — Макрос — Редактор Visual Basic», в открывшемся окне выберем «Insert — Module» и введем:

Function f(ByVal x, ByVal t)

f = 5 * Sin(1 + x) * Exp(2 * t)

End Function

Function mu(ByVal x)

mu = Sin(1 + x)

End Function

Function mu_0(ByVal x)

mu_0 = 2 * Sin(1 + x)

End Function

Function mu_1(ByVal t)

mu_1 = Sin(1) * Exp(2 * t)

End Function

Function mu_2(ByVal t)

mu_2 = Sin(2) * Exp(2 * t)

End Function

Function df(ta, ch, a1, a2, a3, a4)

df = ta * (ch * (a1 - 2 * a2 + a3) + a4)

End Function

Function Ur_Kol_Str(c, a, T1, m, n)

Application. Volatile (True)

Dim u(), y()

ReDim u(n, m), y(n)

h = a / n: tau = T1 / m: h2 = h * h: tau2 = tau ^ 2

c2h2 = (c / h) ^ 2: tau22 = tau2 / 2

For i = 0 To n: u(i, 0) = mu(i * h): Next i

For i = 1 To n - 1

x1 = i * h - h: x2 = i * h: x3 = i * h + h

uxt = df(tau22, c2h2, mu(x1), mu(x2), mu(x3), f(x2, 0))

u(i, 1) = mu(x2) + mu_0(x2) * tau + uxt

Next i

u(0, 1) = mu_1(tau): u(n, 1) = mu_2(tau)

For k = 2 To m

tk = k * tau: tk1 = k * tau - tau: u(0, k) = mu_1(tk): u(n, k) = mu_2(tk)

For i = 1 To n - 1

x2 = i * h

uxt = df(tau2, c2h2, u(i - 1, k - 1), u(i, k - 1), u(i + 1, k - 1), f(x2, tk1))

u(i, k) = 2 * u(i, k - 1) - u(i, k - 2) + uxt

Next i

Next k

For i = 0 To n

y(i) = u(i, m)

Next i

Ur_Kol_Str = y

End Function

Функция Ur_Kol_Str(c, a, T1, m, n) выводит массив-строку из (n + 1)-го элемента. Этот массив равен последнему m-му столбцу матрицы ui,k, которая содержит приближенные значения функции u(x, t).

Чтобы получить результат расчета по этим функциям, выделим диапазон-строку из 11 ячеек, например, B27:L27, введем формулу =Ur_Kol_Str(1;1;1;20;10) и нажмем Ctrl + Schift + Enter. Полученные значения должны совпадать со значениями строки B22:L22.

Из-за громоздкости таблицы расчетов A1:L27 мы её здесь не приводим.

Приведем программу на C++ для решения задачи примера 8.1:

#include <except. h>

#include <iostream. h>

#include <math. h>

long double f(long double x, long double t);

long double mu(long double x);

long double mu_0(long double x);

long double mu_1(long double t);

long double mu_2(long double t);

typedef long double (*PF)(long double, long double);

typedef long double (*PM)(long double);

int Urkolstr(PF f, PM mu, PM mu_0, PM mu_1, PM mu_2, long double c,

long double a, long double T,

long double **u, const int m, const int n);

int main(){

long double **u; long double c, a, T,uxt, uerr; int i, n,m, n1,m1;

PF pf; PM pmu, pmu_0, pmu_1, pmu_2;

cout <<" \n input c = " ; cin >> c;

cout <<" \n input a = " ; cin >> a;

cout <<" \n input T = " ; cin >> T;

cout <<" \n input n = " ; cin >> n; n1 = n+1;

cout <<" \n input m = " ; cin >> m; m1 = m+1;

try {

u = new long double*[n1]; for(i=0;i<=n;i++) u[i]=new long double[m1];

}

catch (xalloc){cout <<"\n Could not allocate\n"; exit(-1);}

pf=f; pmu=mu; pmu_0=mu_0; pmu_1=mu_1; pmu_2=mu_2;

Urkolstr(pf, pmu, pmu_0, pmu_1, pmu_2, c, a, T, u, m, n);

cout << "\n function u: ";

for(i=0;i<=n;i++){cout<<"\n"; uxt=sin(1+i*a/n)*exp(2*T); uerr=fabs((uxt-u[i][m])/uxt);

cout<<" sin(1+x)exp(2)= "<<uxt<<" u["<<i<<"][m]= "<< u[i][m]<<" error= "<<uerr;} cout << "\n Press Key and Enter to continue";

cin >> i; // for pause

for(i=0;i<=n;i++)delete[] u[i];

delete u;

return 0;

}//end main

long double f(long double x, long double t){

long double y;

y = 5*sin(1+x)*exp(2*t);

return y;

}// end f

long double mu(long double x){

long double y;

y = sin(1+x);

return y;

}// end mu

long double mu_0(long double x){

long double y;

y = 2*sin(1+x);

return y;

}// mu_0

long double mu_1(long double t){

long double y;

y = sin(1)*exp(2*t);

return y;

}// mu_1

long double mu_2(long double t){

long double y;

y = sin(2)*exp(2*t);

return y;

}// mu_2

int Urkolstr(PF f, PM mu, PM mu_0, PM mu_1, PM mu_2, long double c,

long double a, long double T,

long double **u, const int m, const int n){

long double h, tau, xi, tk, tau2,c2h2; int i, k;

h=a/n; tau=T/m; tau2=tau*tau; c2h2=c*c/(h*h);