Нахождение приближенного решения краевой задачи гиперболичес

[ewert]

(Оффтоп)

Pphantom в сообщении #1086062 писал(а):

преждевременная оптимизация суть грех.

Не суть, но есть(даже ежели и не есть). Звиняйте за орфографию.

[GAA]

dair, пожалуйста, напишите полностью разностную аппроксимацию задачи. [А не только уравнения и одного из начальных условий. Должны быть учтены все граничные и начальные условия]. Только после этого переходите к попыткам программирования.


[Ваш код в сообщении post1086028.html#p1086028 совсем плохой. Но это станет ясно и Вам самому, как только аккуратно выпишите разностную аппроксимацию всей задачи.]

[dair]

GAA в сообщении #1086081 писал(а):

напишите полностью разностную аппроксимацию задачи

$$w_{h}=\{x_{j}=jh,0..N,hN=1\}, w_{\tau}=\{t_{n}=n\tau, 0...K, K\tau=T\}$
$\frac{y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}}{\tau^2}=1/h((a_{j}-a_{j+1})\frac{y^{n}_{j+1}-y^{n}_{j}}{2h}-(a_{j}-a_{j-1})\frac{y^{n}_{j}-y^{n}_{j-1}}{2h}) + f(x)y^{n}_{j},$$ $j=1...N-1,n=1...K-1$
$y^{n+1}_{0}=0, y^{n+1}_{N}=0, n=0...K-1$
$y^{0}_{j} = 32x_{j}^2(1-x_{j})^3$, $j=1...N-1$
$y_{j}^{1}=y^{0}_j$, $j=1...N-1$

[GAA]

[1. В вашем случае $f(x,t)$ не зависит от $t$, поэтому я убрал в Вашем сообщении неверно набранную зависимость от времени.
2. Вы аппроксимировали второе начальное условие достаточно грубо. Возможно, следует попробовать вариант аппроксимации более высокого порядка по $\tau$.]

Теперь необходимо учесть, что в системе MatLab нумерация индексов начинаются с 1.
Возможно два варианта:

(a) использовать другую нумерацию узлов;
(b) использовать туже нумерацию узлов (по пространственной переменной), но модифицировать разностную аппроксимацию уравнения для узлов 1 и $N-1$ (явно учесть, что в 0 и $N$ значения функции равны нулю).

Для различных шагов по времени обычно изменяют нумерацию: $n = 1$ будет соответствовать начальному моменту времени.

-- Сб 26.12.2015 22:20:23 --

И еще, хоть и не по делу, но, возможно поможет быстрее написать программу.
В вашем случае коэффициент при второй производной не зависит от времени.
Если не менять нумерацию по пространственной координате,
я коэффициенты $a_{j-1/2}$, $j=1…N$ занумеровал бы как $a_j$, $j=1…N$:

$a_j = \cos h(j -1/2) + 2$, $j=1…N.$

Эти коэффициенты вычисляются перед «циклом».
(Перепроверьте нумерацию. Я устал, могу ошибиться.
Конечно, при таком подходе Вам нужно будет убрать 2 в знаменателе аппроксимации уравнения.)

[dair]

Напишу часть кода пока, нумерацию возьму от $1$ до $N+1$

N = 30;
M = 30;
U = zeros(N, M);
tau = 0.01;
h = 0.01;

$y^{n+1}_{0}=0, y^{n+1}_{N}=0, n=0...K-1$

for n=1:M
    U(1,n+1)=0;
    U(N+1,n+1)=0;
end

$y^{0}_{j} = 32x_{j}^2(1-x_{j})^3, j=1...N-1$

for j=2:N
       U(j, 1) = 32*j^2 *(1-j)^3;
end

[GAA]

dair в сообщении #1086117 писал(а):

N = 30;
M = 30;
U = zeros(N, M);
…

for n=1:M
    U(N+1,n+1)=0;
end

Нестыковка: U = zeros(N, M) и U(N+1,n+1)=0.

dair в сообщении #1086117 писал(а):

Напишу часть кода пока, нумерацию возьму от $1$ до $N+1$

N = 30;
h = 0.01;

h = 1/N;

dair в сообщении #1086117 писал(а):

$y^{0}_{j} = 32x_{j}^2(1-x_{j})^3, j=1...N-1$

for j=2:N
       U(j, 1) = 32*j^2 *(1-j)^3;
end

Вам уже писали: не U(j, 1) = 32*j^2 *(1-j)^3, а U(j, 1) = 32*(j*h)^2 *(1-j*h)^3. Но и улучшенный вариант надо проверить. :) [U(j, 1) = 32*(j*h)^2 *(1-j*h)^3 или U(j, 1) = 32*((j-1)*h)^2 *(1-(j-1)*h)^3.]
(Если так будете продолжать и дальше, то тема быстро раздуется, а задание Вы так и не выполните.)

Главное: Вы не привели условие устойчивости Вашей разностной схемы. При любом ли соотношении шагов по пространственной и временной пер. схема будет устойчивой?

[Прим.
1. Обычно первый индекс нумерует шаги по времени, а второй — по пространственной переменной. См., например функцию pdepe.
2. Задание или не очень естественное или что-то подразумевается, но не записано явно.
Шагов по времени может быть очень много. Это может привести к тому, что памяти для хранения всех шагов по времени не хватит. Поэтому обычно сохраняют не все шаги по времени, а часть этих шагов. Сошлюсь в качестве примера на функцию pdepe. При интегрировании шаг по времени выбирается системой автоматически (в отличие от разбиения по пространственной переменной, которое задаёт пользователь). Шаг по времени может быть достаточно малым (для получения требуемой точности). Но сохраняется и возвращается пользователю только значения решения с некоторых временных слоёв. Эта функция MatLab использует не метод конечных разностей, но схема довольно общая.]

[dair]

С учетом Ваших указаний попробую аккуратно записать

N = 30;
M = 30;
U = zeros(M+1, N+1);
tau = 1/M;
h = 1/N;

for n=1:M
    U(n+1,1)=0;
    U(n+1,N+1)=0;
end

for j=2:N
     U(1, j) = 32*(j*h)^2 *(1-j*h)^3;
end
 
for j=2:N
    aj1=cos(h*(j-1/2))+ 2;
    aj2=cos(h*(j+1/2))+ 2;
    fj=sin(pi*j*h);
    for n=2:M
         U(n+1,j) = (tau^2 ./ h)*(aj2*(U(n,j+1)-U(n,j))/h-aj1*(U(n,j)-U(n,j-1))/h)+(tau^2*fj*U(n,j))+2*U(n,j)-U(n-1,j);
    end
end

Еще осталось одно условие, все-таки может оставить оценку, которую я ранее писал, с первым порядком аппроксимации по $\tau$

(Оффтоп)

Если у меня получится выполнить первый пункт, это уже большой упех , так что тема не раздуется надеюсь, а я буду стараться

[GAA]

Так лихо с шагом по времени не получится. Во-первых, по условию интегрировать по времени надо до 5. Во-вторых, устойчивость наложит ограничение на шаг по времени.
(Я третий раз повторять не буду.)

Внешний цикл по времени, а внутренний по пространственной переменной.
Коэффициенты a(j), а не a1(j) и a2(j). В Вашем случае они вычисляются до цикла по времени. В простейшем упражнении это не чувствуется, но если Вам писать что-то серьезное или при защите задания будут внимательно смотреть, то лучше откровенные ляпы не допускать.

-- Вс 27.12.2015 00:26:09 --

dair в сообщении #1086130 писал(а):

С учетом Ваших указаний попробую аккуратно записать

...
for j=2:N
    aj1=cos(h*(j-1/2))+ 2;
    aj2=cos(h*(j+1/2))+ 2;
    fj=sin(pi*j*h);
    for n=2:M
...

for n=3..M+1
M находится по tau и $T=5$.

[dair]

GAA в сообщении #1086134 писал(а):

Так лихо с шагом по времени не получится. Во-первых, по условию интегрировать по времени надо до 5. Во-вторых, устойчивость наложит ограничение на шаг по времени

Надо исследовать наше уравнение $y_{j}^{n+1}=...$ на устойчивость, подставив $y_{j}^{n}=q^{n} e^{ijh\varphi}$ ???
(А еще для уравнения гиперболического типа есть условие Куранта: $\tau \le \min\left(\frac{h}{|\lambda|_{\max}}\right)$)

GAA в сообщении #1086134 писал(а):

M находится по tau и $T=5$.

Найдем условие на $\tau $ и затем $M=T/\tau$
Исправил

   for j=1:N+1 
    a(j)=cos(h*(j-1/2))+ 2;
    f(j)=sin(pi*j*h);
   end
for n=3:M+1
    for j=2:N
         U(n+1,j) = (tau^2 ./ h)*(a(j)*(U(n,j+1)-U(n,j))/h-a(j)*(U(n,j)-U(n,j-1))/h)+(tau^2*f(j)*U(n,j))+2*U(n,j)-U(n-1,j);
    end
end

[GAA]

По поводу устойчивости см., например, главу III (Однородные разностные схемы) книги Самарский А.А. Введение в теорию разностных схем. [Может быть где-то есть проще, чем у Самарского во «Введении…». Я литературой для гиперболических уравнений не интересовался.]

Мелочи:

dair в сообщении #1086137 писал(а):

   for j=1:N+1 
    a(j)=cos(h*(j-1/2))+ 2;
    f(j)=sin(pi*j*h);
   end

В Вашем варианте изменения нумерации узлов, узел j==1 находится на левой границе ($x=0$). Поэтому, f(1) = sin(pi*0), т.е. f(1) = sin(pi*(j-1)*h). Аналогично узел j==N+1 находится на правой границе ($x=1$). Поэтому, f(N+1) = sin(pi), т.е. f(N+1) = sin(pi*(j -1)*h). Следовательно, надо проверить весь цикл вычисления f.

dair в сообщении #1086137 писал(а):

…U(n+1,j) = (tau^2 ./ h)*(a(j)*(U(n,j+1)-U(n,j))/h-a(j)*(U(n,j)-U(n,j-1))/h)+… 

(a(j+1)*(U(n,j+1)-U(n,j))/h-a(j)*(U(n,j)-U(n,j-1))/h)
Перепроверьте индексы.

dair в сообщении #1086137 писал(а):

Найдем условие на $\tau $ и затем $M=T/\tau$

M — целое. Будет немного по-другому. Но это уже совсем мелочи. Сначала — условие устойчивости.

[GAA]

GAA в сообщении #1086134 писал(а):

for n=3..M+1

Это я сбился. Обычно использовал другие обозначения и поздно было; засыпал. В этом месте у Вас как раз было правильно: for n=2:M.

[dair]

GAA в сообщении #1086175 писал(а):

Сначала — условие устойчивости

А можно не проверять условие устойчивости, а подобрать шаг по времени исходя из коэффициентов уравнения?
С устойчивостью пока разбираюсь, но медленно как то идет, поэтому думаю может подобрать

[GAA]

В случае Вашей разностной аппроксимации (постоянные шаги $\tau$, $h$; непрерывный коэффициент при второй производной по пространственной переменной, $k(x)$) условие устойчивости имеет вид

$$\tau < \frac {h}{\sqrt{c_2} }, \qquad c_2 = \max_{x \in [0,1]} k(x).$$

Проверьте по книге Самарского, пожалуйста.

В смысле написания программы, дальше техническая возня, которую надо выполнить аккуратно. Если это семестровое задание, то нужно посмотреть примеры в конспекте лекций или методичке, чтобы сделать то, что от Вас ожидают.

[dair]

Спасибо большое, буду делать!