Волновое уравнение в КР: условие устойчивости и визуализация

[byulent]

Пытаюсь составить решение уравнения колебания струны в конечных разностях и визуализировать его. Собственно, вопрос: как правильно выглядит условие устойчивости для данного уравнения, $\tau < \dfrac{h}{a}$ или $\tau\leqslant\dfrac{h}{a}$?
Если я беру $\tau = \dfrac{h}{2a}$, то график натяжения струны сползает в 0 и дальше не идёт.
Если я беру $\tau = \dfrac{h}{a}$, то график просто перемещается вдоль оси.
Значения на концах: $U^j_1=0, U^j_n=0$.

[Lia]

Просьба более полно сформулировать задачу и расшифровать все свои обозначения, - в противном случае Вы рискуете остаться без ответа.

[Munin]

byulent в сообщении #1070182 писал(а):

график просто перемещается вдоль оси.

Это правильное поведение, так и должно быть. При разных начальных условиях будут и другие варианты поведения.

[byulent]

Munin в сообщении #1070450 писал(а):

byulent в сообщении #1070182 писал(а):

график просто перемещается вдоль оси.

Это правильное поведение, так и должно быть. При разных начальных условиях будут и другие варианты поведения.

А при каком НУ график будет колебаться в привычном нам значении?

[Munin]

Возьмите заранее известное решение "в привычном нам значении", и просто продифференцируйте его по времени. Это и будут начальные условия.

И всё-таки без пояснения ваши обозначения непонятны. Вы полагаете стандартным и общепонятным то, что стандартным и общепонятным не является.

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

[byulent]

Вот смотрите. У меня есть функция $u(x)=px^2+q$, где $p = -4\dfrac{u_{max}}{l^2}, q=4\dfrac{u_{max}}{l}$, $u_{max}$ - максимум функции $u(x)$, $l$ - длина струны. Тогда $U^0_i=u(x_i)$, где $x_i - hi$, $h$ - шаг по $x$, $i=1,2,3\ldots,N$, а, соответственно, $N=\dfrac{l}{h}$.
Считаем, что струна отпущена и $\dfrac{\partial{u}}{\partial{t}}=0$. Тогда $U^1_n=u(x_i)$.
Конечно-разностная схема стандартная:

$$\frac{U^{j-1}_i-2U^j_i+U^{j+1}_i}{\tau^2}=a^2\frac{U^j_{i-1}-2U^j_i+U^j_{i+1}}{h^2},$$

где $\tau$ = шаг по времени.
О граничных условиях и вопросе об условии устойчивости я уже сказал.
Собственно, вид графика с 1 условием устойчивости:

Тут он движется даже, можно сказать, не по оси, а... эммм... сам по себе.
А вот со 2 (правда, тут я заменил 2 на 1.0049):

Упал в 0 и дальше не идёт.
Вот я и спрашиваю, как сделать так, чтобы график не перемещался, а колебался нормально, то есть постепенно бы превращался в своё отражение в минусе, а затем опять уходил в плюс.

[Munin]

Выглядит как что-то неправильно запрограммированное. Потому что картинка в любом случае должна быть симметричной по горизонтали.

[byulent]

Munin в сообщении #1077763 писал(а):

Выглядит как что-то неправильно запрограммированное. Потому что картинка в любом случае должна быть симметричной по горизонтали.

Что насчёт программирования...

Используется синтаксис C++

class PhysObject
{
public:
    PhysObject(int len);
    ~PhysObject();
    double* getBeg();
    double* getNext(double *a, int n);
    int getN();
    int getU();
    double getDx();
    double getDt();
private:
    static const double dx = 0.01;
    static const double A = 1;
    static const int U = 1;
    double dt;
    int n;
    double *Ubeg;
    double *Uprev;
    double *Unext;
    double* UAllocate (int n);
};

PhysObject::PhysObject(int len)
{
    n = len/dx;
    dt = dx/(1.0049*A);
    Uprev = UAllocate(n);
    Ubeg = new double [n+1];
    for (int i=0;i<=n;i++) Ubeg[i] = Uprev[i];
    //std::cout<<std::endl;
    Unext = new double [n+1];
}

double* PhysObject::UAllocate(int n)
{
    double *a = new double[n+1];
    for (int i=0; i<=n; i++){
        double L= n*dx;
        double p=-4*U/(L*L);
        double q=4*U/L;
        double x=i*dx;
        a[i]=p*x*x+q*x;
        //a[i]=5*sin(x);
        //std::cout<<a[i] << " ";
    }
    //std::cout<<std::endl;
    return a;
}

double* PhysObject::getNext(double *a, int n)
{
    for (int i=1; i<n-1; i++){
        double k = pow(A,2)*pow(dt,2)/pow(dx,2);
        Unext[i]=k*(a[i+1]-2*a[i]+a[i-1])+2*a[i]-Uprev[i];
        Uprev[i]=a[i];
        //std::cout<<Unext[i] << " ";
    }
    Unext[0] = 0;
    Unext[n-1] = 0;
    //std::cout<<std::endl;
    return Unext;
}

int PhysObject::getN()
{
    return n;
}

int PhysObject::getU()
{
    return U;
}

double PhysObject::getDt()
{
    return dt;
}

double PhysObject::getDx()
{
    return dx;
}

PhysObject::~PhysObject()
{
    delete Ubeg;
    delete Uprev;
    delete Unext;
}

double *PhysObject::getBeg()
{
    return Uprev;
}

Что здесь неправильно?

[Munin]

Сдаюсь, копаться в программе мне точно лень.

[Pphantom]

byulent в сообщении #1077827 писал(а):

Что здесь неправильно?

А то его знает, копаться в таком коде тяжело. Но я могу подтвердить предположение Munin - проблема именно в реализации. Нормальная программа, реализующая Ваше описание, выдает ожидаемые результаты в обоих случаях. Собственно в описании есть только одна ошибка:

byulent в сообщении #1077754 писал(а):

меня есть функция $u(x)=px^2+q$,

Должно быть $u(x) = p x^2+q x$ (в программе, кстати, это сделано правильно).

[byulent]

Нагуглил разницу между бегущими и стоячими волнами. Похоже, моя проблемы скрывалась в этом: я хотел получить стоячую волну, а получилась бегущая. Сейчас буду думать, как смоделировать стоячую.

[Pphantom]

byulent в сообщении #1078241 писал(а):

Сейчас буду думать, как смоделировать стоячую.

Уже подумали. Еще раз повторяю - у Вас все верно в описании алгоритма, ищите ошибку непосредственно в программе.

[byulent]

А что здесь может быть не так? Я делал всё по алгоритму. Если вы говорите, что вам якобы сложно разбирать код, то тут собственно вся модель находится в 3 методах: конструкторе, UAllocate и getNext.

[Pphantom]

byulent в сообщении #1078283 писал(а):

А что здесь может быть не так? Я делал всё по алгоритму./quote] Кто ж его знает... Делать-то Вы делали, но явно в процессе где-то ошиблись.

byulent в сообщении #1078283 писал(а):

Если вы говорите, что вам якобы сложно разбирать код, то тут собственно вся модель находится в 3 методах: конструкторе, UAllocate и getNext.

И что? Код малочитаем (об использовании ООП при численном моделировании я помолчу - цензурных слов нет), и лезть в него никому не хочется.

А так... к сообщению прицеплена анимация - результат счета по Вашему алгоритму для случая $\tau = \frac{h}{2a}$. У Вас должно получаться нечто аналогичное.

[byulent]

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

byulent в сообщении #1078283 писал(а):

А что здесь может быть не так? Я делал всё по алгоритму./quote] Кто ж его знает... Делать-то Вы делали, но явно в процессе где-то ошиблись.

byulent в сообщении #1078283 писал(а):

Если вы говорите, что вам якобы сложно разбирать код, то тут собственно вся модель находится в 3 методах: конструкторе, UAllocate и getNext.

И что? Код малочитаем (об использовании ООП при численном моделировании я помолчу - цензурных слов нет), и лезть в него никому не хочется.

А так... к сообщению прицеплена анимация - результат счета по Вашему алгоритму для случая $\tau = \frac{h}{2a}$. У Вас должно получаться нечто аналогичное.

У меня получается аналогичное, но график движется жутко медленно и доходит только до 0. Кроме того, иногда он может ВНЕЗАПНО подскочить в $+\infty$. При повторной прогонке с этими же значениями всё происходит нормально.