2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2, 3, 4, 5 ... 7  След.
 
 численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 15:28 
Интеграл табличный, нужно для проверки правильности ввода.
$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} {\exp(-ax^2-by^2-cx-dy)} dxdy$. Файл: http://zalil.ru/34871895/5dea9e66.52bcc620/int2.nb. Выдается куча ошибок
$
H= 7.2\cdot10^5;
\Theta0 = 0.018;
\Delta f = 3.5\cdot10^8;
P (t) = NSum[
  NIntegrate [
   Exp[-((5.55(y^2 + x^2))/(
      H^2\Theta0^2)) - \pi \Delta f^2(t - x^2/(
        cH) - y^2/(cH))^2], {x, -\infty, 
    \infty}, {y, -\infty, \infty}, MaxRecursion \to40, 
   AccuracyGoal \to60], {k, -255, 255}]$

 
 
 
 Re: численноt решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 15:36 
Аватара пользователя
Неохота скачивать файл. Напишите здесь выражение.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 16:33 
salang
Так, начнём с внутреннего выражения (числ. интегрирования). У вас там не задано $\[c\]$ и $\[t\]$, о чём система прямо и сообщает. Вы хотя бы читаете её сообщения? Далее. Вы пытаетесь суммировать по k, причём k у вас в выражении отсутствует.
P.S.И присоединюсь к просьбе ИСН, выкладывайте код сюда под соотв. тегом.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 16:53 
вот полный файл (предыдущий я упростил): http://zalil.ru/34872046/6e965a53.52bcd5c0/6.nb.
А как выложить в сообщение? Там кучу исправлений надо внести.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 17:10 
1)У вас не объявлена переменная "da".
2)Использовать в индексе переменной такой же индекс, как и по тому, что вы суммируете не хорошо (если это изначально так не задумано) - у вас это $\[{\sigma _k}\]$.
3)Неудачная конструкция сумма(интеграл).
4)Впрочем это не главная проблема. Если отдельно брать ваши интегралы (для разных k), система говорит, что они крайне близки к 0 в любые из указанных вами моментов времени. Тут нужно аккуратно разобраться. Во первых, явно можно сделать кучу замен и сильно упростить выражение, и потом уже " в читабельном виде" разбираться.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 17:14 
Аватара пользователя
salang, вносите все эти исправления и вставляйте код прямо в сообщение под тегом [code].

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 17:55 
п.1 и 2 исправил- http://zalil.ru/34872154/3c37fa27.52bce560/6.nb
п.3 не понял
п.4 все выражения, содержащие тригонометрию уже заменены на экспоненты, не совсем понятно, что еще можно упростить
Код:
c = 3*10^8;
H = 7.2*10^5;
\[Tau] = 47*10^-6;
Fd = 7200;
Tp = 50*10^-6;
\Delta f = 3.5*10^8;
v = 7500;
Q = 18000;
L = 100;
Subscript[\[Sigma], 2] = 0.003;
Subscript[\[Sigma], 1] = 2;
\[CapitalTheta]0 = 0.018;
m = 256;
\[Lambda] = 0.022;
n = 3;
lm = 25;
d = 1.2;
\[Lambda] = 0.022;
n = 3;
t = List[ -10^-8, -10^-9, 0, 10^-16, 10^-13, 10^-9, 10^-8, 20*10^-8];

P (t) = NSum[
  NIntegrate [
   Exp[x^2*(\[Pi]*Fd^2*\[Tau]^2*(m - Abs[k])^2 )/H^2 - (
     0.00005*n^2*\[Pi]^2*d^2*y^2 )/(\Theta0^2*\[Lambda]^2 *H^2) - (
     100*lm^2*x^2*Subscript[\[Sigma], 1]^2)/(
     H^2*L^2*\[Pi]^2*Subscript[\[Sigma], 2]^2) - (
     100*lm^2*y^2*Subscript[\[Sigma], 1]^2)/(
     H^2*L^2*\[Pi]^2*Subscript[\[Sigma], 2]^2) - (5.55*(y^2 + x^2))/(
     H^2*\Theta0^2) - \[Pi] * (\Delta f^2 *(t - k*Tp - x^2/(
           c*H) - y^2/(c*H))^2 +
        2*Q*(v*k*Tp)/(H*d) + \[Tau]^2*((v*k*Tp)/(H*d))^2 )], {x, -Infinity,
    Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
   AccuracyGoal -> 60], {k, -255, 255}]

ListLinePlot[%, AxesLabel -> {"t", "P(t)"}, PlotStyle -> PointSize[0.01]]

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 20:22 
Аватара пользователя
Ох, тут скорее не упрощать, а просто исправлять надо.
1. С какой целью используется NSum? Это одна из ошибок, замените на простое Sum.
2. Если интеграл численно не хочет считаться методом по умолчанию, то можно ему явным образом задать метод, какой выбирайте из справки.
3. "P(t)=" - такая запись в математике не годится, уберите ее, тем более что у Вас это вообще не функция t? Впрочем это как раз на результат особо не влияет.
4. Вам точно нужны суммы интегралов по каждому t? Я просто уточняю на всякий случай.

В общем вот такой код вполне себе будет работать
Код:
Sum[NIntegrate[
  Exp[x^2*(\[Pi]*Fd^2*\[Tau]^2*(m - Abs[k])^2)/
      H^2 - (0.00005*n^2*\[Pi]^2*d^2*
       y^2)/(\[CapitalTheta]0^2*\[Lambda]^2*H^2) - (100*lm^2*x^2*
       s1^2)/(H^2*L^2*\[Pi]^2*s2^2) - (100*lm^2*y^2*s1^2)/(H^2*
       L^2*\[Pi]^2*
       s2^2) - (5.55*(y^2 +
         x^2))/(H^2*\[CapitalTheta]0^2) - \[Pi]*(df^2*(t - k*Tp -
           x^2/(c*H) - y^2/(c*H))^2 +
       2*Q*(v*k*Tp)/(H*
           d) + \[Tau]^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
   Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
  AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255, 255}]

Я тут некоторые символам заменил названия, мне так больше нравится.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 20:47 
В идеальном случае хорошо бы получить непрерывную кривую от времени, если это вообще возможно для моих данных. Дискретное время я задал просто для упрощения расчета. P на самом деле является непрерывной функцией времени и там есть зависимость от него, остающаяся и после интегрирования.
C последним кодом тоже не получилось: http://zalil.ru/34872426/5b108a0.52bd1058/6.nb

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 21:21 
Аватара пользователя
Ну вы даете... Я же написал, что поменял символы, привел только последнюю часть кода, разумеется без изменения первой части она работать не будет, переопределите те три константы, которые я обозвал по своему и все будет работать.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение26.12.2013, 22:13 
не выходит каменный цветок: http://zalil.ru/34872573/672b068d.52bd23e0/6.nb
Код:
Sum[NIntegrate[
  Exp[x^2*(\[Pi]*Fd^2*\[Tau]^2*(m - Abs[k])^2)/
      H^2 - (0.00005*n^2*\[Pi]^2*d^2*
       y^2)/(\[CapitalTheta]0^2*\[Lambda]^2*H^2) - (100*lm^2*x^2*
       s1^2)/(H^2*L^2*\[Pi]^2*s2^2) - (100*lm^2*y^2*s1^2)/(H^2*
       L^2*\[Pi]^2*
       s2^2) - (5.55*(y^2 +
         x^2))/(H^2*\[CapitalTheta]0^2) - \[Pi]*(f^2*(t - k*Tp -
           x^2/(c*H) - y^2/(c*H))^2 +
       2*Q*(v*k*Tp)/(H*
           d) + \[Tau]^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
   Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
  AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255, 255}]

511 NIntegrate[
  Exp[(x^2 (\[Pi] Fd^2 \[Tau]^2 (m - Abs[k])^2))/H^2 - (
    0.00005 n^2 \[Pi]^2 d^2 y^2)/(\[CapitalTheta]0^2 \[Lambda]^2 H^2) \
- (100 lm^2 x^2 s1^2)/(H^2 L^2 \[Pi]^2 s2^2) - (100 lm^2 y^2 s1^2)/(
    H^2 L^2 \[Pi]^2 s2^2) - (5.55 (y^2 + x^2))/(
    H^2 \[CapitalTheta]0^2) - \[Pi] (f^2 (t - k Tp - x^2/(c H) - y^2/(
          c H))^2 + (2 Q (v k Tp))/(
       H d) + \[Tau]^2 ((v k Tp)/(
         H d))^2)], {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \
\[Infinity]}, MaxRecursion -> 40, AccuracyGoal -> 60,
  Method -> "AdaptiveMonteCarlo"]

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение28.12.2013, 13:46 
Аватара пользователя
salang, сосредоточьтесь и выполните последовательно весь код из Вашего файла.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение28.12.2013, 22:47 
я сделал это сразу же после Вашего сообщения. Но что-то не получается: http://zalil.ru/34875236/12247ab.52bfcb90/6.nb. Может просто выложить Ваш рабочий вариант?

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение28.12.2013, 22:53 
Аватара пользователя
Вы лучше выложите тот вариант, который пробовали последним. И пожалуйста, приводите код прямо в сообщении, это ведь несложно.

 
 
 
 Re: численное решение интеграла от экспоненты в Mathematica
Сообщение29.12.2013, 09:18 
запросто
Код:
In[49]:= c = 3*10^8;
H = 7.2*10^5;
\[Tau] = 47*10^-6;
Fd = 7200;
Tp = 50*10^-6;
f = 3.5*10^8;
v = 7500;
Q = 18000;
L = 100;
s2 = 0.003;
s1 = 2;
\[CapitalTheta]0 = 0.018;
m = 256;
\[Lambda] = 0.022;
n = 3;
lm = 25;
d = 1.2;
\[Lambda] = 0.022;
n = 3;


In[69]:= P (t) = Sum[
  NIntegrate [
   Exp[x^2*(\[Pi]*Fd^2*\[Tau]^2*(m - Abs[k])^2 )/H^2 - (
     0.00005*n^2*\[Pi]^2*d^2*y^2 )/(\[CapitalTheta]0^2*\[Lambda]^2 *H^2) - (
     100*lm^2*x^2*s1^2)/(H^2*L^2*\[Pi]^2*s2^2) - (100*lm^2*y^2*s1^2)/(
     H^2*L^2*\[Pi]^2*s2^2) - (5.55*(y^2 + x^2))/(
     H^2*\[CapitalTheta]0^2) - \[Pi] * (f^2 *(t - k*Tp - x^2/(c*H) - y^2/(
           c*H))^2 +
        2*Q*(v*k*Tp)/(H*d) + \[Tau]^2*((v*k*Tp)/(H*d))^2 )], {x, -Infinity,
    Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
   AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255, 255}]


During evaluation of In[69]:= Set::write: Tag Times in P {-(1/100000000),-(1/1000000000),0,1/1000000000,1/100000000,1/10000000,1/1000000,1/5000} is Protected. >>

Out[69]= {1.59372*10^-12, 343561., 626592., 964461., 59001.6, 0., 0., 676401.}

In[70]:= ListLinePlot[%, AxesLabel -> {"t", "P(t)"}, PlotStyle -> PointSize[0.01]]

Out[70]= \!\(\*
GraphicsBox[{{}, {{}, {},
{RGBColor[0.24720000000000017`, 0.24, 0.6], PointSize[0.01],
      LineBox[{{1., 1.5937237677846204`*^-12}, {2., 343560.79215414973`}, {3.,
        626591.6960827865}, {4., 964461.1900936544}, {5., 59001.578511843}, {
       6., 0.}, {7., 0.}, {8., 676401.3207484136}}]}}, {}},
AspectRatio->0.6180339887498948,
Axes->True,
AxesLabel->{
FormBox["\"t\"", TraditionalForm],
FormBox["\"P(t)\"", TraditionalForm]},
AxesOrigin->{0, 0},
Method->{},
PlotRange->{{0, 8.}, {0, 964461.1900936544}},
PlotRangeClipping->True,
PlotRangePadding->{{0.16, 0.16}, {19289.22380187309, 19289.22380187309}}]\)

In[45]:= Sum[NIntegrate[
  Exp[x^2*(\[Pi]*Fd^2*\[Tau]^2*(m - Abs[k])^2)/
      H^2 - (0.00005*n^2*\[Pi]^2*d^2*y^2)/(\[CapitalTheta]0^2*\[Lambda]^2*
       H^2) - (100*lm^2*x^2*s1^2)/(H^2*L^2*\[Pi]^2*s2^2) - (100*lm^2*y^2*
       s1^2)/(H^2*L^2*\[Pi]^2*
       s2^2) - (5.55*(y^2 +
         x^2))/(H^2*\[CapitalTheta]0^2) - \[Pi]*(f^2*(t - k*Tp - x^2/(c*H) -
           y^2/(c*H))^2 +
       2*Q*(v*k*Tp)/(H*d) + \[Tau]^2*((v*k*Tp)/(H*d))^2)], {x, -Infinity,
   Infinity}, {y, -Infinity, Infinity}, MaxRecursion -> 40,
  AccuracyGoal -> 60, Method -> "AdaptiveMonteCarlo"], {k, -255, 255}]

Out[45]= {1.5983*10^-12, 335063., 836068., 936968., 90404., 0., 0., 702430.}

In[46]:= ListLinePlot[%, AxesLabel -> {"t", "P(t)"}, PlotStyle -> PointSize[0.01]]

Out[46]= \!\(\*
GraphicsBox[{{}, {{}, {},
{RGBColor[0.24720000000000017`, 0.24, 0.6], PointSize[0.01],
      LineBox[{{1., 1.598299392712959*^-12}, {2., 335062.66072612564`}, {3.,
       836067.9196629806}, {4., 936968.342828458}, {5., 90403.98299123366}, {
       6., 0.}, {7., 0.}, {8., 702429.7247402478}}]}}, {}},
AspectRatio->0.6180339887498948,
Axes->True,
AxesLabel->{
FormBox["\"t\"", TraditionalForm],
FormBox["\"P(t)\"", TraditionalForm]},
AxesOrigin->{0, 0},
Method->{},
PlotRange->{{0, 8.}, {0, 936968.342828458}},
PlotRangeClipping->True,
PlotRangePadding->{{0.16, 0.16}, {18739.36685656916, 18739.36685656916}}]\)

t = List[ -10^-8, -10^-9, 0, 10^-9, 10^-8, 10^-7, 10^-6, 20*10^-5];

Как сделать плавную кривую на графике? Массив значений t удалил, но результат все равно не очень.
Как сделать, чтобы по x на графике показывались не номера отсчетов, а абсолютное время?

 
 
 [ Сообщений: 93 ]  На страницу 1, 2, 3, 4, 5 ... 7  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group