Метод проекционной регуляризации

[vortex2]

Доброго времени суток, уважаемые форумчане!

Пусть дана обратная задача нахождения функции $h(t)$ по имеющимся значениям $f(t)$

$$ \frac{\partial u_1(x,t)}{\partial t}=a_1^2\frac{\partial^{2}v_1(x,t)}{\partial x^2},\ \ 0<x<x_0,\ \ t>0, $$

$$ \frac{\partial u_2(x,t)}{\partial t}=a_2^2\frac{\partial^{2}u_2(x,t)}{\partial x^2},\ x_0<x<1,\ \ t>0,\ \ a_1,\ a_2>0, $$

$$ u_1(x,0)=0, \ \ x\in[0,x_0];\qquad u_2(x,0)=0,\ \ x\in \left( \right. x_0,1 \left. \right], $$

$$ u_1(0,t)=h(t),\ \ t\geq 0, $$

$$ u_2(1,t)=0,\ \ \ t\geq0, $$

$$ u_1(x_0,t)=u_2(x_0,t),\ \ t\geq 0, $$

$$ a_1\frac{\partial u_1(x,t)}{\partial x}=a_2\frac{\partial u_2(x,t)}{\partial x},\ \ t\geq 0, $$

$$u_1(x_0,t)=f(t),\quad t\geq 0.$$

И также имеется алгоритм численного решения методом проекционной регуляризации с исходными данными

$$ a_1,~a_2,~x_0,~T,~f_{\delta}(t),~\delta. $$

1. Применяем к $f_{\delta}(t)$ преобразование Фурье, находим $\hat{f_{\delta}}(\tau)$.

2. Нахождение $ \hat{h}_{\delta}^{\overline{\alpha}(\hat{f}_{\delta}(\tau),\delta)}=\hat{h}_{\delta}(\tau) $ методом проекционной регуляризации

$$ \hat{h}^{\alpha}_\delta(\tau)= \begin{cases} \frac{\sh\mu_0\biggl(\frac{(1-x_0)\sqrt{\tau}}{a_2}+\frac{x_0\sqrt{\tau}}{a_1}\biggr)}{\sh\mu_0(1-x_0)\frac{\sqrt{\tau}} {a_2}}\cdot\hat{f}_\delta(\tau), \mu_0=(1+i)/\sqrt{2}, & \tau \leq \alpha,\\ 0, & \tau > \alpha. \end{cases} $$

Значение параметра регуляризации $\alpha$ получим из принципа невязки

$$ \int_{\alpha}^{\infty}|\hat{f}_{\delta}(\tau)|^2d\tau=9\delta^2. $$

Обозначим полученное значение параметра регуляризации через $\overline{\alpha}(\hat{f}_{\delta}(\tau),\delta)$.

3. Применяя обратное преобразование Фурье $F^{-1}$, получаем приближенное решение $h_\delta(t)$

$$ h_\delta(t)= \begin{cases} \operatorname{Re}[F^{-1}[\hat{h}_\delta(\tau)]], & t \in [0,T],\\ 0, & t \notin [0,T]. \end{cases} $$

В качестве модельного примера была выбрана функция $ h_0(t)= t\sin(\pi t), & t \in [0,1]. $

Попробовал реализовать, вот что получилось (не судите строго)

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

public alglib.complex OperatorT(double r, double a1, double a2, double x0, double tau)
{ 
    double arg2 = (1 - x0) * Math.Sqrt(tau) / a2 / Math.Sqrt(2);
    double arg1 = arg2 + x0 * Math.Sqrt(tau) / a1 / Math.Sqrt(2);

    alglib.complex T = new alglib.complex();
    if (tau == 0) {
        T.x =  (x0 * a2 + a1 + x0 * a1) / (a1 - x0 * a1); // находится с помощью предела
        T.y =  (x0 * a2 + a1 + x0 * a1) / (a1 - x0 * a1);
  
    }
    else if (tau > r) { T.x = 0; T.y = 0; } // применение параметра регуляризации
    else
    {
        T.x = (Math.Sinh(arg1) * Math.Cos(arg1) * Math.Sinh(arg2) * Math.Cos(arg2) 
            + Math.Cosh(arg1) * Math.Sin(arg1) * Math.Cosh(arg2) * Math.Sin(arg2)) 
            / (Math.Sinh(arg2) * Math.Cos(arg2) * Math.Sinh(arg2) * Math.Cos(arg2) 
            + Math.Cosh(arg2) * Math.Sin(arg2) * Math.Cosh(arg2) * Math.Sin(arg2));
        T.y = (Math.Cosh(arg1) * Math.Sin(arg1) * Math.Sinh(arg2) * Math.Cos(arg2) 
            - Math.Sinh(arg1) * Math.Cos(arg1) * Math.Cosh(arg2) * Math.Sin(arg2)) 
            / (Math.Sinh(arg2) * Math.Cos(arg2) * Math.Sinh(arg2) * Math.Cos(arg2) 
            + Math.Cosh(arg2) * Math.Sin(arg2) * Math.Cosh(arg2) * Math.Sin(arg2));
    }
    return T;
}

private void ButtonS_Click(object sender, RoutedEventArgs e)
{
    
    int n = 512; // число узлов разбиения временного отрезка
    
    alglib.xparams _params = new alglib.xparams(9);
    
    alglib.complex [] a1 = new alglib.complex[n];
    alglib.complex [] a2 = new alglib.complex[n];
    alglib.complex[] a3 = new alglib.complex[n];
    
    double [] r1 = new double[n];
    double [] r2 = new double[n];
    double[] r3 = new double[n];
    
    double t = 1; // длина временного отрезка
    
    double h = (double)t / (n - 1); //шаг разбиения
    
    double[] xS = new double[n];
    double[] Fs = new double[n];
    
    double delta = 0.01; // уровень погрешности в начальных данных
    
    Random r = new Random();
    
    for (int i = 0; i <= n - 1; i++)
    {
        xS[i] = i * h;
        r1[i] = xS[i] * Math.Sin(Math.PI * xS[i]); // функция из модельного примера
    }
    double x0 = 0.5, aa1 = 1, aa2 = 1; // ввод начальных данных
    
    Par(r1, x0, ref r2); // нахождение r2 из решения прямой задачи в точке x0
        // r2 - значения функции f(t), используемой в качестве начальных данных для обратной задачи 
    for (int i = 0; i <= n - 1; i++)
    { // внесение шума в начальные данные
        double temp = -1 + 2 * r.NextDouble();
        double dt = delta * temp / Math.Sqrt(t);
        a1[i].x = r2[i] + dt;
        a1[i].y = 0;
    }
    alglib.fft.fftc1d(ref a1, n, _params); // прямое преобразование Фурье

    double td = 200; // длина отрезка по тау
    double th = (double)td / (n - 1); // шаг разбиения
    double R = 6;    // параметр регуляризации    

    for (int i = 0; i <= n - 1; i++)
    { // оператор в методе проекионной регуляризации
        alglib.complex T2 = new alglib.complex();
        T2 = OperatorT(R, aa1, aa2, x0, i * th);

       a2[i].x = (a1[i].x * T2.x - a1[i].y * T2.y); 
       a2[i].y = (a1[i].x * T2.y + a1[i].y * T2.x);
       
    }
    
    alglib.fft.fftc1dinv(ref a2, n, _params); // обратное преобразование Фурье

    for (int i = 0; i <= n - 1; i++)
    {
        r3[i] = a2[i].x; // выделение Re
    }

    // отображение графиков
    var X = xS.AsXDataSource();
    var Y = r3.AsYDataSource();
    var FS = r1.AsYDataSource();

    CompositeDataSource compositeDataSource = X.Join(Y);
    CompositeDataSource compositeDataSourceF = X.Join(FS);       
    plotter.AddLineGraph(compositeDataSource, Colors.Green, 2, "приближенное");
    plotter.AddLineGraph(compositeDataSourceF, Colors.Red, 2, "точное");


}
 

Вот что получается



Меняя параметр регуляризации R, а также длину $\tau$ ничего лучше получить не удается. Само собой такое решение никуда не годится. Возможно я что-то делаю принципиально неправильно. Помогите, пожалуйста, разобраться.

Может есть примеры реализации подобной задачи, например, на MATLAB?

[thething]

vortex2 в сообщении #1488678 писал(а):

Само собой такое решение никуда не годится.

Почему такой вывод? Есть какие-то оценки погрешности, в рамки которых Вы не попадаете?

vortex2 в сообщении #1488678 писал(а):

Возможно я что-то делаю принципиально неправильно.

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

[vortex2]

thething в сообщении #1488689 писал(а):

Почему такой вывод? Есть какие-то оценки погрешности, в рамки которых Вы не попадаете?

Должно получиться наподобии такого (здесь разве что модельная задача немного другая):



В моём случае даже без оценок погрешности видно, что с решением что-то не так.

thething в сообщении #1488689 писал(а):

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

В данном случае у меня вопросы по реализации. Оператор $R\hat{f}(\tau)=\sh\mu_0\biggl(\frac{(1-x_0)\sqrt{\tau}}{a_2}+\frac{x_0\sqrt{\tau}}{a_1}\biggr)\cdot \sh^{-1}\mu_0(1-x_0)\frac{\sqrt{\tau}}{a_2}\cdot\hat{f}(\tau)$ найден верно.

[thething]

thething в сообщении #1488689 писал(а):

Должно получиться наподобии такого (здесь разве что модельная задача немного другая):

Правильно я понимаю, с другой модельной задачей Ваша программа сработала хорошо?

А насчёт оценок -- это Вы зря. В таких задачах погрешности решения могут во много раз превосходить погрешности исходных данных. Визуальное сравнение тут -- не показатель.

vortex2 в сообщении #1488694 писал(а):

В данном случае у меня вопросы по реализации

Тут ничем помочь не могу, ибо си-шарп не знаю от слова совсем.

Аналогичные задачи с программными реализациями надо, по-видимому искать в диссертациях, в открытом доступе такие программы вряд ли есть.

[vortex2]

thething в сообщении #1488700 писал(а):

Правильно я понимаю, с другой модельной задачей Ваша программа сработала хорошо?

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

[Pphantom]

Если я правильно понял код (C# я тоже не знаю), у вас там вносится дополнительный шум, полная амплитуда которого раза в три больше, чем у "сигнала". Вы уверены, что в такой ситуации результат должен хоть сколько-нибудь соответствовать ожидаемому? Что будет, если зашумление полностью убрать?

[vortex2]

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

Если я правильно понял код (C# я тоже не знаю), у вас там вносится дополнительный шум, полная амплитуда которого раза в три больше, чем у "сигнала". Вы уверены, что в такой ситуации результат должен хоть сколько-нибудь соответствовать ожидаемому? Что будет, если зашумление полностью убрать?

Если полностью убрать шум, то почти ничего не меняется



А вот что будет, если вместо R=6 задать R=60 (без шума и с шумом соответственно).

[sup]

vortex2 в сообщении #1488694 писал(а):

В данном случае у меня вопросы по реализации. Оператор

$$ R\hat{f}(\tau)=\sh\mu_0\biggl(\frac{(1-x_0)\sqrt{\tau}}{a_2}+\frac{x_0\sqrt{\tau}}{a_1}\biggr)\cdot \sh^{-1}\mu_0(1-x_0)\frac{\sqrt{\tau}}{a_2}\cdot\hat{f}(\tau) $$

найден верно.

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

$$ \sinh (X) \cosh (Y) + \frac{a_2}{a_1}\sinh (Y) \cosh (X) $$

Мне кажется, что Вы потеряли множитель $\frac{a_2}{a_1}$. После чего числитель свернулся в простой синус.

-- Сб окт 24, 2020 08:50:48 --

Хотя, там же условия склейки другие.
Виноват.

-- Сб окт 24, 2020 09:24:25 --

Еще один вопрос. Вы используете $FFT$. Это же означает, что вместо решения начально-краевой задачи, Вы ищете решение, периодическое по $t$. Для того, чтобы снизить искажения, возможно, стоит продолжить $f(t)$ нулем на бОльший интервал и решать там.

Для проверки собственно схемы вычислений можно попробовать в качестве $f(t)$ брать чистые экспоненты (комплексные). Для них все промежуточные вычисления можно явно проконтролировать. Только лучше брать показатель "несоизмеримый" с интервалом $(0,T)$.

[vortex2]

sup в сообщении #1488786 писал(а):

Еще один вопрос. Вы используете $FFT$. Это же означает, что вместо решения начально-краевой задачи, Вы ищете решение, периодическое по $t$. Для того, чтобы снизить искажения, возможно, стоит продолжить $f(t)$ нулем на бОльший интервал и решать там.

Спасибо. Да, действительно, после увеличения интервала получилось что-то похожее на правду.



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

[sup]

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

$$ \begin{align} &u_t - u_{xx} = 0, &(t, x) \in (0,T) \times (0,1) \\ &u(0,x) = 0, &x \in (0,1),\\ &u(t,0) = 0, &t \in (0,T),\\ &u_x(t,0) = f(t), &t \in (0,T), \end{align} $$

Требуется найти $v(t) = u(t, 1)$. Т.е. вместо того, чтобы задавать краевые условия (по $x$) на обеих стенках, задаем данные Коши на одной стенке.
Легко показать, что эта задача эквивалентна уравнению Вольтерры первого рода сверточного типа

$$ \int \limits_0^t K(t - s)v(s) \, ds = f(t). $$

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

Если уж Вы хотите использовать $FFT$, то лучше это сделать для данного интегрального уравнения.
Ядро $K(t)$ можно предъявить в виде некого "неудобного" ряда. А можно явно указать преобразование Лапласа $\hat K(t)$.