Matlab символическая интеграция с реальным дает сложный ответ - PullRequest
0 голосов
/ 09 декабря 2018

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

close all;
clear all;
clc;

syms t real;
syms thetak real;
syms sik real;
syms tbar real;
syms sjk real;

expr = exp(-thetak*((t-sik)^2 + (t-sjk)^2));
Bijk_raw = int(expr,t,0,1);
Bijk = simplify(collect(expand(Bijk_raw)));
fprintf('Bijk is as follows...\n');
pretty(Bijk);

Ответы [ 2 ]

0 голосов
/ 09 декабря 2018

Ответ, который вы получите (если у вас есть версия Matlab, похожая на мою) и которую я воспроизвожу здесь:

  /               /                     2 \ 
  |  1/2   1/2    |   thetak (sik - sjk)  | 
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2          1/2                 \ 
 |    | 2    (-thetak)    (sik i + sjk i) | 
 | erf| --------------------------------- | i - 
 \    \                 2                 / 

    /  1/2          1/2                       \   \ \ 
    | 2    (-thetak)    (sik i + sjk i - 2 i) |   | | 
 erf| --------------------------------------- | i | | / 
    \                    2                    /   / / 
             1/2 
 (4 (-thetak)   )

создает впечатление, что у вас есть комплексное число i везде.

Но на самом деле это ложное впечатление из-за (-thetak) ^ (1/2).

Действительно, получение квадратного корня из отрицательного числа сгенерирует «i», которое в свою очередь «убьет»другие "я", которые "соприкасаются" с ним.Эта отмена будет происходить в разных местах из-за того, что (-thetak) ^ (1/2) можно найти:

1) внутри выражений erf и

2) как общеезнаменатель (последняя строка).

Убедитесь, что правило i ^ 2 = -1 применяется везде, не оставляя шансов на выживание любого "i" ...

Наконец дайте (я установил thetak= s ^ 2 с s> 0):

  /               /                      \ 
  |  1/2   1/2    |   s^2 (sik - sjk)^2   |   
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2                   \ 
 |    | 2    s    (sik  + sjk ) | 
 | erf| ----------------------- |  - 
 \    \           2            / 

    /  1/2                          \   \ \ 
    | 2    s   (sik  + sjk  - 2 )   |   | | 
 erf| ----------------------------- |   | |   /   (4 s)
    \             2                 /   / / 

Редактировать: Вы могли избежать интеграции.Идея состоит в том, чтобы преобразовать квадратик в $ exp (-thetak * ((t-sik) ^ 2 + (t-sjk) ^ 2)) $ в так называемую «каноническую форму», которая в вашем случае: $ exp(-thetak * (((tA) ^ 2 + B)) / C); $, где $ A, B, C $ можно выразить как функции sik и sjk (например, $ A = (sik + sjk) / 2$);таким образом, установив $ T = tA $, вы вернетесь к классическому интегралу Гаусса с формулой:

$$ \ frac {2} / {\ sqrt {\ pi}} \ int_a ^ b exp(-t ^ 2} dt) (erf (b) - erf (a)) $$

0 голосов
/ 09 декабря 2018

То, что вы получаете, имеет (с некоторыми постоянными коэффициентами) форму

 1i * (c * erf(1i * a) - erf(c * 1i * (a - 2)))

Состоит из двух слагаемых формы

- 1i * erf(1i * x)

, которая также известна как мнимая функция ошибок erfi().Оказывается, что

erfi(x) = - 1i * erf(1i * x) = 2/sqrt(pi) * integral(@(t)exp(t.^2),0,x)

Так что для реальных значений x ваше выражение на самом деле является реальным, и это имеет место, если thetak >= 0 и sik и sjk являются действительными числами.

Интеграл, с которого вы начинаете, может быть уменьшен до интеграла exp(-t^2) (с использованием некоторого аффинного преобразования), который, как известно, не имеет "замкнутой формы", но обычно записывается как

erf(x) = 2/sqrt(pi) * integral(@(t)exp(-t.^2),0,x)

Я настоятельно рекомендую прочитать статью в Википедии о функции .

Кроме того, я рекомендую использовать более удобную CAS для новичков, чем символьная панель инструментов MATLAB.Один бесплатный CAS с открытым исходным кодом, который я хотел бы порекомендовать, это Maxima .

(Это все написано в нотации MATLAB из-за отсутствия LaTeX здесь на SO.)

...