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

//1.

for (i=0; i<=n; i++) u[i][0]=mu(i*h); // k=0

u[0][1]=mu_1(tau); u[n][1]=mu_2(tau);

for (i=1; i<n; i++){xi = i*h; u[i][1]=mu(xi)+mu_0(xi)*tau+

tau2*(c2h2*(mu(xi-h)-2*mu(xi)+mu(xi+h))+f(xi,0));} // k=1

//2. k=2,...,m-1: u[i][k]:

for(k=2; k<=m;k++){ tk=k*tau; u[0][k]=mu_1(tk); u[n][k]=mu_2(tk);

for (i=1; i<n; i++){xi=i*h;

u[i][k]=2*u[i][k-1]-u[i][k-2]+

tau2*(c2h2*(u[i-1][k-1]-2*u[i][k-1]+u[i+1][k-1])+f(xi, tk-tau));}

} // end k

return 0;

}// Urkolstr

Приведем листинг работы программы.

Как видим, результаты расчета совпадают до трех знаков после запятой со значениями точного решения для t = T = 1.

8.2.2.  Уравнение колебаний струны. Неявная схема

Построим неявную схему для уравнения колебаний струны (8.11):

(8.25)

Для устойчивости схемы (8.25) параметры c, h, τ, σ должны удовлетворять условию [9]

.

Если ¼ ≤ σ ≤ 1/2, то схема (8.25) безусловно устойчива.

Шаблон схемы (8.25) изображен на рис. 8.3.

Рис. 8.3

Значения решения на нулевом и первом слоях вычисляют по формулам (8.22), (8.23). На каждом k-ом слое (k = 2, 3, …, M) решают методом прогонки систему уравнений относительно :

(8.26)

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

0. Построить сеточную область, выбирая шаги h, τ и параметр σ так, чтобы выполнялось условие устойчивости

.

Если ¼ ≤ σ ≤ 1/2, то можно выбирать h, τ произвольно.

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

(8.27)

(8.28)

2. Значения ui,k+1 для каждого k > 0 находим методом прогонки, последовательными вычислениями в несколько этапов.

2.1. Вычислим правые части (8.26):

(8.29)

2.2. Вычислим прогоночные коэффициенты:

(8.30)

(8.31)

(8.32)

2.3. Вычислим решение ui,k+1:

(8.33)

(8.34)

Пример 8.2. Решить задачу из примера 8.1, используя неявную схему.

Решение на C++. Приведем текст программы, реализующей алгоритм (8.27) — (8.34):

#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 sgm,

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

int main(){

long double **u; long double c, a, T,sgm, 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;

cout <<" \n input sgm (1/4 <= sgm <= 1/2) = " ; cin >> sgm;

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, sgm, 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(2T)= "<<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 sgm,

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

long double *fi,*ai,*bi, h,tau, xi, tk, tau2,c2h2,h2c2t2s, h2c2s, sgm1,d; int i, k,n1;

h=a/n; tau=T/m; tau2=tau*tau/2; c2h2=c*c/(h*h);h2c2t2s=1/(c2h2*tau*tau*sgm);

n1=n+1; d=2+h2c2t2s; h2c2s=1/(c2h2*sgm); sgm1=(1-2*sgm)/sgm;

//1.

for (i=0; i<=n; i++) u[i][0]=mu(i*h); //(8.27) k=0

u[0][1]=mu_1(tau); u[n][1]=mu_2(tau);

for (i=1; i<n; i++){xi = i*h; u[i][1]=mu(xi)+mu_0(xi)*tau+

tau2*(c2h2*(mu(xi-h)-2*mu(xi)+mu(xi+h))+f(xi,0));} // (8.28) k=1

//2. k=1,2,...,m-1: u[i][k+1]:

try {fi = new long double [n1];ai = new long double [n1];

bi = new long double [n1];}

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

for(k=1; k<m;k++){ tk=k*tau;

//2.1

for (i=1; i<n; i++){xi=i*h; fi[i]=h2c2t2s*(-2*u[i][k]+u[i][k-1])-

sgm1*(u[i-1][k]-2*u[i][k]+u[i+1][k])-

(u[i-1][k-1]-2*u[i][k-1]+u[i+1][k-1])-h2c2s*f(xi, tk);}

//2.2

ai[1]=1/d; bi[1]=(fi[1]-mu_1(tk+tau))/(-d);

for(i=2; i<n-1; i++){ai[i]=-1/(ai[i-1]-d);

bi[i]=(fi[i]-bi[i-1])/(ai[i-1]-d); }

bi[n-1]=(fi[n-1]-mu_2(tk+tau)-bi[n-2])/(ai[n-2]-d);

//2.3

u[0][k+1]=mu_1(tk+tau); u[n][k+1]=mu_2(tk+tau);

u[n-1][k+1]=bi[n-1];

for(i=n-2; i>0; i--)u[i][k+1]=ai[i]*u[i+1][k+1]+bi[i];

} // end k

delete[] bi;

delete[] ai;

delete[] fi;

return 0;

}// Urkolstr

Приведем листинг работы программы:

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

Отметим преимущества неявной схемы перед явной схемой:

В явной схеме надо выбирать шаги h и τ так, чтобы выполнялось условие устойчивости (условие Куранта) cτ < h. В неявной схеме можно выбирать шаги h и τ произвольно, если параметр σ принадлежит отрезку
[1/4; 1/2].

8.2.3.  Разностный метод для уравнения колебаний мембраны

Рассмотрим задачу для уравнения колебаний однородной прямоугольной мембраны с закрепленными краями

(8.37)

(8.38)

(8.39)

Введем сеточную область

Обозначим . Заменяя производные разностными формулами, для уравнения (8.37) получим разностное уравнение с порядком аппроксимации :

(8.40)

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

Решение на нулевом слое k = 0 получим из начального условия (8.38):

(8.41)

Чтобы найти решение на первом слое (k = 1) используем разложение в ряд Тейлора по переменной t в окрестности t = 0:

(8.42)

Учитывая начальные и краевые условия получим приближенные значения решения на первом слое:

Для вычисления значений на следующих слоях (k > 0) имеем из (8.40):

(8.43)

Приведем условие устойчивости явной схемы (8.41) — (8.43) [9]:

(8.44)