Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Numerical Convolution Problem, different results by Mathematica V3 and V5.0

  • To: mathgroup at smc.vnet.net
  • Subject: [mg73436] Numerical Convolution Problem, different results by Mathematica V3 and V5.0
  • From: "Zhao, Liang" <ZhaoL at MedImmune.com>
  • Date: Fri, 16 Feb 2007 01:04:12 -0500 (EST)
  • References: <64F125B1D61F974FA9CF8362865ACE7E0F52E3@MD1EV002.medimmune.com>

Hello,
I guess I am running into a numerical problem when I am trying to 
perform Convolutions on InterpolatingFunctions .
Firstly, I get the "exact" answer by performing the convolution with 
closed form functions and assign values to integral limits as shown as 
follows.

test=DSolve[{x'[t]==y[t], y'[t]==x[t],x[0]==1, 
y[0]==2}, {x[t], y[t]}, t]

ff1[t_]=test[[1]][[1]][[2]]

ff2[t_]=test[[1]][[2]][[2]]

convolve[f_,g_,t_]:=Integrate[f[u]*g[t-u],{u,0,t}]

ff3[t_]=convolve[ff1,ff2,t]

N[ff3[10]]

ff4[t_]=N[convolve[ff3,ff1,t]]

N[ff4[10]]

ff5[t_]:=convolve[ff4,ff3, t]

N[ff5[10]]

The above code gives:

495595.
3.5972 10^6

6.53123 10^7

However, when I calculate the same thing from the numerical route, it 
gives different result.
< /FONT>< /FONT>

test = NDSolve[{x'[t] == y[t],

       y'[t] == x[t], x[0] == 1, y[0] == 2}, {x[t], y[t]}, 
{t, 0, 10}];

f1[t_] = test[[1]][[1]][[2]];

f2[t_] = test[[1]][[2]][[2]];

convolve[f_, g_, t_?NumericQ] := N[NIntegrate[f[u]*g[t - u], {u, 0, 
t}]];

f3[t_?NumericQ] := convolve[f1, f2, t];

f3[10]

f4[t_?NumericQ] := convolve[f3, f1, t];

f4[10]

f5[t_?NumericQ] := convolve[f4, f3, t];

f5[10]

It gives

495595.

4.79075 10^6

5.13256 10^17

Even more interestingly, when I sent the same code over to friend to run 
it with Mathematica V3.0, the numerical solution rendered by the above 
code is fine. I need your help to unpuzzle the myth!


Liang Zhao



  • Prev by Date: Re: Quick integral.
  • Next by Date: Re: Showing that a hypergeometric expression is 0?
  • Previous by thread: Dll problem
  • Next by thread: How to apply a function to a function multiple times