使用NDSolve的输出作为NIntegrate的积分限制


4

我一直在研究ODE,对于函数y[t],由于它的复杂性,必须以数字方式求解。然后我需要使用NDSolve的结果,这是我存放在变量 $s$ ,以下列方式:

Plot[ 
  NIntegrate[ 
    f[r],{r,b,Evaluate[Sqrt[{(y[t])^2+b^2}/.s]]} 
    ], 
  {t,0,5} 
  ] 

返回以下错误:

NIntegrate::nlim: r = {20.025} is not a valid limit of integration. >> 
NIntegrate::nlim: r = {20.025} is not a valid limit of integration. >> 
NIntegrate::nlim: r = {19.6024} is not a valid limit of integration. >> 
General::stop: Further output of NIntegrate::nlim will be suppressed during this calculation. >> 

看来这个问题是由于 $text{Evaluate}$ 返回括号中的数字......我尝试过 $text{N}$ 这个函数,但它有同样的问题。也许我在做别的事情。谁能帮我?

编辑:我将张贴在该意见提出的完整代码:

Module[{b = 1., k = 1., Vt0 = 2., L0 = 20.}, 
s = NDSolve[ 
  {y''[t] - (y'[t])^2 (b^2/((b^2 + (y[t])^2) y[t]) - y[t]/(
   b^2 (Exp[(y[t])^2/b^2] - 1))) - (z'[t])^2 (-1/(
   y[t] Exp[(y[t])^2/b^2]) (Exp[(y[t])^2/b^2] - 
   1) (b^2 + (y[t])^2)) + k*(b^2 y[t])/(b^2 + (y[t])^2) == 0, 
  y[0] == L0, y'[0] == 0, 
  z''[t] + 2*y'[t]*z'[t]*(y[t]/(b^2 + (y[t])^2)) == 0, z'[0] == Vt0, 
  z[0] == 0}, 
  {z, y}, {t, 0, 5} 
  ]; 
f[r_] := Sqrt[E^(r^2/b^2)/(-E + E^(r^2/b^2)) - 1]; 

Plot[ 
  NIntegrate[f[r], {r, b, Evaluate[Sqrt[(y[t])^2 + b^2] /. s]}], 
  {t, 0, 5} 
  ] 
] 
  0

你应该把信息放在'f [],s,r' .. !! 25 8月. 132013-08-25 12:26:39

  0

@纳西尔ups。是\ sqrt是Sqrt []。在使用'code'命令之前,我复制为Latex,并且保留了这个工件。 25 8月. 132013-08-25 12:29:54

  0

@Blackbird只需编辑以包含完整的代码。函数$ \ text {z [t]} $与错误无关。 25 8月. 132013-08-25 12:33:24

  0

'M'说你要指定一个插值函数,它需要一个常量。你需要研究这个。可能你的方法是错误的。'NIntegrate'的上限不能是函数返回函数。 25 8月. 132013-08-25 12:44:13

  0

@Blackbird根据你的建议,我找出了我的方法出了什么问题。非常感谢您的宝贵时间。 25 8月. 132013-08-25 14:03:08

4

我不知道,如果这个情节使得没有意义。但是你的复杂值从你的NIntegrate中显示出来。任何方式,可能你可以看看这个,看看它看起来不错。

让我们分散一点。一旦它工作正常,那么你可以创建一个模块并将其折叠成一个很好的长表达式。

b = 1.; 
k = 1.; 
Vt0 = 2.; 
L0 = 20.; 

sol = NDSolve[{y''[ 
     t] - (y'[t])^2 (b^2/((b^2 + (y[t])^2) y[t]) - 
     y[t]/(b^2 (Exp[(y[t])^2/b^2] - 1))) - (z'[ 
      t])^2 (-1/(y[t] Exp[(y[t])^2/b^2]) (Exp[(y[t])^2/b^2] - 
      1) (b^2 + (y[t])^2)) + k*(b^2 y[t])/(b^2 + (y[t])^2) == 0, 
    y[0] == L0, y'[0] == 0, 
    z''[t] + 2*y'[t]*z'[t]*(y[t]/(b^2 + (y[t])^2)) == 0, z'[0] == Vt0, 
    z[0] == 0}, {z, y}, {t, 0, 5}]; 
ysol = y /. [email protected]; 

确保解决方案的下一步

Plot[ysol[t], {t, 0, 5}] 

Mathematica graphics

f[r_?NumericQ] := Sqrt[E^(r^2/b^2)/(-E + E^(r^2/b^2)) - 1]; 

检查

f[r] 

Mathematica graphics

012之前看起来确定

使数据。我用Abs,因为你有复杂的数字。

data = Table[{t,Abs @ NIntegrate[f[r], {r, b, Sqrt[ysol[t]^2 + b^2]}]}, {t, 0,5, .05}]; 

检查几点:

data[[1 ;; 3]] 

Mathematica graphics

情节

ListLinePlot[data, Joined -> True, Mesh -> True, MeshStyle -> Red, 
PlotRange -> All, AxesLabel -> {"t", "NIntegrate"}] 

Mathematica graphics

另一种方式做上面是先计算一个向量的上限,确保它没问题,然后再用下一个进行积分。发现相同的结果。

upperLimit = Table[Sqrt[ysol[t]^2 + b^2], {t, 0, 5, .01}]; 
ListLinePlot[upperLimit, Mesh -> True, MeshStyle -> Red] 

Mathematica graphics

data = Table[{i, Abs @ NIntegrate[f[r], {r, b, upperLimit[[i]]}]}, {i, 
    Length[upperLimit]}]; 

ListLinePlot[data, Joined -> True, Mesh -> True, MeshStyle -> Red, 
PlotRange -> All, AxesLabel -> {"i", "NIntegrate"}] 

Mathematica graphics

我没有看到什么有趣的在这里。让t步更小,就可以得到直线。

  0

我相信你的方法是正确的。 我已经知道我的问题了。我只需要将上限表达式等于一个变量,例如$ \ text {x} $。然后在NIntegrate中使用$ \ text {x [[1]]} $作为上限。但是当然你的方法非常有趣,我预见我会经常检查它。非常感谢你 25 8月. 132013-08-25 14:01:29