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

[dair]

Помогите написать код в Matlab-e или подскажите с чего начать. Задание:

Цитата:

1. Напишите программу в MatLab для поиска приближенного решения краевой задачи для одномерного уравнения гиперболического типа:

$$ u_{tt} = \left((\cos x +2)u_x)_x - (\sin \pi x)u, \qquad u=u(x,t), \qquad 0 \le x \le1, \qquad 0 \le t \le 5, $$

$$ u(0,t) = u(1,t) = 0, \qquad u(x,0) = 32x^2(1-x)^3, \qquad u_t(x,0) = 0.$$

Стройте график решения для каждого момента времени $t$.
2. Постройте график функции $E(t)$:

$$E(t) = \int\limits_0^t (u_t^2+u_x^2) dx.$$

Для численного вычисления интеграла используйте метод прямоугольников.
3. Постройте в координатах $(x,t)$ трехмерную поверхность $u(x,t)$. Примечание. Графики необходимо выводить в одно графическое окно, но в разные его части. Не допускается строить различные графики в одной плоскости. Графики должны быть оформлены должным образом: оси координат должны быть подписаны, должна быть нанесена сетка, у каждого графики должен быть заголовок.

Нашел пример реализации явной разностной схемы для уравнениях

$u_y = u_{xx} + 8x^2y^4$, $u(x,0) = \sin(\pi x)$, $-u_x(0,y) + 0.5 u(0,y) = 0$, $u_x(1,y) + 0.5u(1,y) = 0.$

function result=yavn(M)
pause;
N = 3*M^2
h=1/M;
t=1/N;
U=zeros(M+1,N+1);
for X=1:M+1,
   U(X,1)=sin(pi*(X-1)*h);
end
for J=1:N,
   y=(J-1)*t;    
   U(1,J+1)=2*t/h^2*(U(2,J)-U(1,J)*(1+h/2-h^2/(2*t)));
   U(M,J+1)=2*t/h^2*(U(M+1,J)+U(M,J)*(1+h/2-h^2/(2*t)))+t*8*h^2*y^4;
   for X=2:M,
      x=(X-1)*h;
      U(X,J+1)=U(X,J)+t/h^2*(U(X-1,J)-2*U(X,J)+U(X+1,J))+t*(8*x^2*y^4);
   end
end
mesh(0:(1/N):1,0:(1/M):1,U),grid;
result=U;

[GAA]

Найденный Вами пример... не в тему. [В нем рассматривается уравнение параболического типа.]

Приведите, пожалуйста, разностную схему для Вашей задачи. Если Вы не знаете, как составить разностную схему, посмотрите книги по численным методам. [Можно поискать ссылки в темах этого форума. Например в прилепленной теме. В частности, можно посмотреть: Самарский А.А., Гулин А.В. Численные методы; Березин И.С., Жидков Н.П. Методы вычислений, Т2.]

[dair]

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

Найденный Вами пример... не в тему. [В нем рассматривается уравнение параболического типа.]

Приведите, пожалуйста, разностную схему для Вашей задачи. Если Вы не знаете, как составить разностную схему, посмотрите книги по численным методам. [Можно поискать ссылки в темах этого форума. Например в прилепленной теме. В частности, можно посмотреть: Самарский А.А., Гулин А.В. Численные методы; Березин И.С., Жидков Н.П. Методы вычислений, Т2.]

Прочитал Гулина ЧМ про "Разностные схемы для уравнения колебаний"
Используем сетку
$w_{h}=\{x_{i}=ih,0..N,hN=1\}, w_{\tau}=\{t_{n}=n\tau, 0...K, K\tau=T\}$
Надо сделать замену $y^{n}_{tt,j}=((y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}))/\tau^2$ и $y^{n}_{xx,j}=((y^{n}_{j+1}-2y^{n}_{j}+y^{n}_{j-1}))/h^2$ для $j=1...N-1,n=1...K-1$ и использовать начальное условие $y^{0}_{i}=u_{0}(x_{i}), i=1...N-1$, получим разностную схему, а вот как это проделать на моем примере не совсем ясно, подкажите

[Pphantom]

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

а вот как это проделать на моем примере не совсем ясно

Давайте начнем с такого вопроса:

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

Надо сделать замену $y^{n}_{tt,j}=((y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}))/\tau^2$

Как получено это выражение?

[dair]

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

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

а вот как это проделать на моем примере не совсем ясно

Давайте начнем с такого вопроса:

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

Надо сделать замену $y^{n}_{tt,j}=((y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}))/\tau^2$

Как получено это выражение?

Используем пятиточечный шаблон(крест) в котором используется три слоя $n-1,n,n+1$ и запишем
$y_{t}=\frac{y^{n+1}_{j}-y^{n}_{j}}{\tau}$
$y_{\bar{t}}=\frac{y^{n}_{j}-y^{n-1}_{j}}{\tau}$
$y_{\hat{t}}=\frac{y^{n+1}_{j}-y^{n-1}_{j}}{2\tau}$
$y^{n}_{\bar{t}t,j}=\frac{y_{t}-y_{\bar{t}}}{\tau}=\frac{((y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}))}{\tau^2}$
аналогично для $y^{n}_{\bar{x}x,j}$, а как заменить именно вот это $((\cos{x}+2)u_{x})_{x}$?

-- 26.12.2015, 13:34 --

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

а как заменить именно вот это $((\cos{x}+2)u_{x})_{x}$?

Напишу свой вариант: Пусть $a_{i}=\cos{x_{i}+2}$ тогда запишем разностное отношение $(au_{\bar{x}})_{x}=1/h(a_{i+1}\frac{y^{n}_{i+1}-y^{n}_{i}}{h}-a_{i}\frac{y^{n}_{i}-y^{n}_{i-1}}{h}) $

[ewert]

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

запишем разностное отношение $(au_{\bar{x}})_{x}=1/h(a_{i+1}\frac{y^{n}_{i+1}-y^{n}_{i}}{h}-a_{i}\frac{y^{n}_{i}-y^{n}_{i-1}}{h})
$

Только $a_{i\pm\frac12}$

[dair]

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

$y^{n}_{\bar{t}t,j}=\frac{y_{t}-y_{\bar{t}}}{\tau}=\frac{((y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j}))}{\tau^2}$

ewert в сообщении #1085939 писал(а):

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

запишем разностное отношение$ (au_{\bar{x}})_{x}=1/h(a_{i+1}\frac{y^{n}_{i+1}-y^{n}_{i}}{h}-a_{i}\frac{y^{n}_{i}-y^{n}_{i-1}}{h}) $

Только $a_{i\pm\frac12}$

Тогда разностная схема запишется так!?
$\frac{(y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j})}{\tau^2}=1/h(a_{i+1/2}\frac{y^{n}_{i}-y^{n}_{i}}{h}-a_{i-1/2}\frac{y^{n}_{i}-y^{n}_{i-1}}{h}) + f(x,j)y^{n}_{j} $
$u(x,0) = 32x_{i}^2(1-x_{i})^3 $

[Pphantom]

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

Тогда разностная схема запишется так!?
$\frac{(y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j})}{\tau^2}=1/h(a_{i+1/2}\frac{y^{n}_{i}-y^{n}_{i}}{h}-a_{i-1/2}\frac{y^{n}_{i}-y^{n}_{i-1}}{h}) + f(x,j)y^{n}_{j}$

Да, именно. Единственная добавка: $a_{i+1/2}=(a_i+a_{i+1})/2$ (в том смысле, что интерполировать стоит уже косинус, а не считать эту функцию в $x_{i+1/2}$).

[dair]

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

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

Тогда разностная схема запишется так!?
$\frac{(y^{n+1}_{j}-2y^{n}_{j}+y^{n-1}_{j})}{\tau^2}=1/h(a_{i+1/2}\frac{y^{n}_{i}-y^{n}_{i}}{h}-a_{i-1/2}\frac{y^{n}_{i}-y^{n}_{i-1}}{h}) + f(x,j)y^{n}_{j}$

Да, именно. Единственная добавка: $a_{i+1/2}=(a_i+a_{i+1})/2$ (в том смысле, что интерполировать стоит уже косинус, а не считать эту функцию в $x_{i+1/2}$).

Я правильно понял, что теперь надо выразить
$ y^{n+1}_{j}=\tau^2/h((a_j+a_{j+1})\frac{y^{n}_{i}-y^{n}_{i}}{2h}-(a_j+a_{j-1})\frac{y^{n}_{i}-y^{n}_{i-1}}{2h}) + \tau^2f(x,j)y^{n}_{j}+2y^{n}_{j}-y^{n-1}_{j}$
и начинать с $y^{0}_{j}= 32x_{j}^2(1-x_{j})^3$ до $j=N-1$ и $n=K-1 $

-- 26.12.2015, 17:25 --

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

$a_{i+1/2}=(a_i+a_{i+1})/2$

что то не совсем понятно, откуда вы это получили пожалуйста поподробней объясните

[Pphantom]

Кстати, один индекс $i+1$ Вы потеряли (в первом слагаемом, при $y$). И индексы где-то $i$, где-то $j$, хотя должна быть какая-то одна буква.

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

Я правильно понял, что теперь надо выразить

Да - с точностью до поправок выше.

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

что то не совсем понятно, откуда вы это получили

А как Вы хотели бы вычислять что-то на полуцелой сетке?

[ewert]

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

Единственная добавка: $a_{i+1/2}=(a_i+a_{i+1})/2$ (в том смысле, что интерполировать стоит уже косинус, а не считать эту функцию в $x_{i+1/2}$).

Нет никакого смысла интерполировать косинус (хоть и приведёт это к практически тем же результатам), коль скоро по условию он считается тупо явно.

-- Сб дек 26, 2015 19:33:20 --

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

что то не совсем понятно, откуда вы это получили

В результате тупой линейной интерполяции. Которая вполне сохраняет второй порядок точности схемы в целом; но которая вполне и не нужна.

[dair]

Попробовал написать программу, что то ошибки повалились, может не так написал

N = 30;
M = 30;
U = zeros(N, M);
tau = 0.01;
h = 0.01;
 
 
% численный метод
for k=1:N-1
 U(k, 1) = 32*k^2 *(1-k)^3;
end
 
 
for n=1:14
    for i=1:M-1
         U(i, n + 1) = (tau^2 ./ h)*((cos(i)+cos(i+1)+4)*(U(i+1,n)-U(i,n))/2*h-(cos(i)+cos(i-1)+4)*(U(i,n)-U(i-1,n))/2*h)+(tau^2*sin(pi*i)*U(i,n))+2*U(i,n)-U(i,n-1);
    end
end
 

[Pphantom]

ewert в сообщении #1086025 писал(а):

Нет никакого смысла интерполировать косинус (хоть и приведёт это к практически тем же результатам), коль скоро по условию он считается тупо явно.

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

-- 26.12.2015, 19:29 --

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

что то ошибки повалились, может не так написал

У Вас аргументом того самого пресловутого косинуса (да и синуса) служит номер узла. Какой в этом смысл? Не путайте координату узла и его номер.

[ewert]

(Оффтоп)

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

Не стоит учить человека плохому.

И в то же время стОит учить товарища грамотному, в смысле Оккаму.

[Pphantom]

(Оффтоп)

ewert в сообщении #1086056 писал(а):

И в то же время стОит учить товарища грамотному, в смысле Оккаму.

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