Is it possible to numerically solve a differential equation by FFT?

First order ODE,

it's of the form

y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),

y(0)=0,

where f(.) may be a nonlinear function,

How to use FFT to do it?

Thanks a lot! 

Mike wrote:

> First order ODE,

> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),

> y(0)=0,
> 
> where f(.) may be a nonlinear function,

> How to use FFT to do it?

FFT works well for two point boundary value problems
which are otherwise difficult to solve.

For a first order initial value problem there are
many ways to do it.

-- glen

On 19 May, 05:02, "Mike" <meathea...@gmail.com> wrote:
> First order ODE,
>
> it's of the form
>
> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),
>
> y(0)=0,
>
> where f(.) may be a nonlinear function,
>
> How to use FFT to do it?
>
> Thanks a lot!


I am not familiar with FFT transforms in this respect, but before
embarking on more complicated procedures, I would probably try to
numerically integrate it by approximating y(t) by a simple step
function and then iterating the equation. So if you rewrite your
equation as

(y(t)-y(t-dt))/dt =c1*y(t-dt) +c2 +c3*exp(f(y(t-dt)) ,  i.e.

y(t)= y(t-dt) + dt*[c1*y(t-dt) +c2 +c3*exp(f(y(t-dt))]

you can calculate y(t) from your previous value y(t-dt) (for instance
for y(0)=0, you would get y(dt)=dt*(c2+c3*exp(f(0))). Plugging then
your new value as the old value into the right hand side gives you
your next value, and so on.
It may not be the most efficient way of integrating the equation, but
it is probably the best way to keep control over your results as there
is only one parameter required, namely dt. Just by making the latter
smaller, it is easy to check if your solution has sufficiently
converged or not.

Thomas

"Thomas Smid" <thomas.smid@gmail.com> wrote in message 
news:1179757633.099586.27890@y18g2000prd.googlegroups.com...
> On 19 May, 05:02, "Mike" <meathea...@gmail.com> wrote:
>> First order ODE,
>>
>> it's of the form
>>
>> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),
>>
>> y(0)=0,
>>
>> where f(.) may be a nonlinear function,
>>
>> How to use FFT to do it?
>>
>> Thanks a lot!
>
>
> I am not familiar with FFT transforms in this respect, but before
> embarking on more complicated procedures, I would probably try to
> numerically integrate it by approximating y(t) by a simple step
> function and then iterating the equation. So if you rewrite your
> equation as
>
> (y(t)-y(t-dt))/dt =c1*y(t-dt) +c2 +c3*exp(f(y(t-dt)) ,  i.e.
>
> y(t)= y(t-dt) + dt*[c1*y(t-dt) +c2 +c3*exp(f(y(t-dt))]
>
> you can calculate y(t) from your previous value y(t-dt) (for instance
> for y(0)=0, you would get y(dt)=dt*(c2+c3*exp(f(0))). Plugging then
> your new value as the old value into the right hand side gives you
> your next value, and so on.
> It may not be the most efficient way of integrating the equation, but
> it is probably the best way to keep control over your results as there
> is only one parameter required, namely dt. Just by making the latter
> smaller, it is easy to check if your solution has sufficiently
> converged or not.
>
> Thomas
>

Thanks a lot Thomas. Are you suggesting a home-made approach to solve the 
ODE?

If it works it would be great! My headache is that I got something out from 
Matlab's ODE45, and ODE 23, and they all give warning of dividing by zero 
and the results have NANs... I am not sure why... and want another approach 
which is guaranteed to be numerically stable and correct so I could compare 
the results ... 

"Mike" <meatheadIV@gmail.com> writes:

> 
> "Thomas Smid" <thomas.smid@gmail.com> wrote in message 
> news:1179757633.099586.27890@y18g2000prd.googlegroups.com...
> > On 19 May, 05:02, "Mike" <meathea...@gmail.com> wrote:
> >> First order ODE,
> >>
> >> it's of the form
> >>
> >> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),
> >>
> >> y(0)=0,
> >>
> >> where f(.) may be a nonlinear function,
> >>
> >> How to use FFT to do it?

> If it works it would be great! My headache is that I got something out from 
> 
> Matlab's ODE45, and ODE 23, and they all give warning of dividing by zero 
> and the results have NANs... I am not sure why... and want another approach 
> 
> which is guaranteed to be numerically stable and correct so I could compare 
> 
> the results ... 

Without knowing the precise function f, I would suspect that either the
solution is blowing up (y -> infinity) at a finite value of t, namely
t = int_0^infty 1/(c1 y + c2 + c3 exp(f(y))) dy, or it's hitting a value
of y where f is undefined.

For example, try y' = y + 1 + exp(y), y(0) = 0; the solution goes to 
infinity at approximately t= .5871365611.
-- 
Robert Israel              israel@math.MyUniversitysInitials.ca
Department of Mathematics        http://www.math.ubc.ca/~israel 
University of British Columbia            Vancouver, BC, Canada

On 22 May, 03:48, "Mike" <meathea...@gmail.com> wrote:
> "Thomas Smid" <thomas.s...@gmail.com> wrote in message
>
> news:1179757633.099586.27890@y18g2000prd.googlegroups.com...
>
>
>
> > On 19 May, 05:02, "Mike" <meathea...@gmail.com> wrote:
> >> First order ODE,
>
> >> it's of the form
>
> >> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),
>
> >> y(0)=0,
>
> >> where f(.) may be a nonlinear function,
>
> >> How to use FFT to do it?
>
> >> Thanks a lot!
>
> > I am not familiar with FFT transforms in this respect, but before
> > embarking on more complicated procedures, I would probably try to
> > numerically integrate it by approximating y(t) by a simple step
> > function and then iterating the equation. So if you rewrite your
> > equation as
>
> > (y(t)-y(t-dt))/dt =c1*y(t-dt) +c2 +c3*exp(f(y(t-dt)) ,  i.e.
>
> > y(t)= y(t-dt) + dt*[c1*y(t-dt) +c2 +c3*exp(f(y(t-dt))]
>
> > you can calculate y(t) from your previous value y(t-dt) (for instance
> > for y(0)=0, you would get y(dt)=dt*(c2+c3*exp(f(0))). Plugging then
> > your new value as the old value into the right hand side gives you
> > your next value, and so on.
> > It may not be the most efficient way of integrating the equation, but
> > it is probably the best way to keep control over your results as there
> > is only one parameter required, namely dt. Just by making the latter
> > smaller, it is easy to check if your solution has sufficiently
> > converged or not.
>
> > Thomas
>
> Thanks a lot Thomas. Are you suggesting a home-made approach to solve the
> ODE?
>
> If it works it would be great! My headache is that I got something out from
> Matlab's ODE45, and ODE 23, and they all give warning of dividing by zero
> and the results have NANs... I am not sure why... and want another approach
> which is guaranteed to be numerically stable and correct so I could compare
> the results ...

Hi Mike,

Yes, I would definitely use a separate integration routine (where I
know exactly how it works) in order to be sure that the results your
are getting are correct (or if you can't get a solution at all). I
wrote a little C++ program according to my suggestion above, which
should do this:


void main() {

  ofstream outfile ("out.txt");

  double y,tmax,t,dt,dtprint;

  tmax=10.0;
  dt=0.001;
  dtprint=0.1;

  t=0;
  y=0;                 /* y(t=0) */

  outfile <<setprecision (4) <<t<<' '<<y<<'\n';


  int nmax=tmax/dt;
  int nprint=tmax/dtprint;


  for (int n=1; n<=nmax; n=n+1)  {
  t=t+dt;
  y=y+dt*(YOUR EXPRESSION FOR y'(t) GOES HERE, WITH y REPLACING y(t));
  if (n%nprint==0)    { outfile <<setprecision (4) <<t<<'
'<<y<<'\n'; }
  }


system("PAUSE");
return 0;

}


I tested the program for the equation y'(t)=c*y(t) (for which the
solution is y(t)=exp(c*t)) and it gave the correct result to within an
accuracy of about 1% with the present stepwidth dt=0.001 (for more
accurate results or for functions that vary more strongly, you'll have
to choose a smaller value for dt).
Note also that you may have to choose an initial value of y different
from zero in order to avoid problems with undefined functions (use a
small value like 1E-10 for instance in these cases).

Thomas

In article <f2tlkl$pe7$1@news.Stanford.EDU>,
 "Mike" <meatheadIV@gmail.com> writes:
 >
 >"Thomas Smid" <thomas.smid@gmail.com> wrote in message 
 >news:1179757633.099586.27890@y18g2000prd.googlegroups.com...
 >> On 19 May, 05:02, "Mike" <meathea...@gmail.com> wrote:
 >>> First order ODE,
 >>>
 >>> it's of the form
 >>>
 >>> y'(t)=c1*y(t) + c2 + c3*exp(f(y(t)),
 >>>
 >>> y(0)=0,
 >>>
 >>> where f(.) may be a nonlinear function,
 >>>
 >>> How to use FFT to do it?
 >>>
 >>> Thanks a lot!
 >>
 >>
 >> I am not familiar with FFT transforms in this respect, but before
 >> embarking on more complicated procedures, I would probably try to
 >> numerically integrate it by approximating y(t) by a simple step
 >> function and then iterating the equation. So if you rewrite your
 >> equation as
 >>
 >> (y(t)-y(t-dt))/dt =c1*y(t-dt) +c2 +c3*exp(f(y(t-dt)) ,  i.e.
 >>
 >> y(t)= y(t-dt) + dt*[c1*y(t-dt) +c2 +c3*exp(f(y(t-dt))]
 >>
 >> you can calculate y(t) from your previous value y(t-dt) (for instance
 >> for y(0)=0, you would get y(dt)=dt*(c2+c3*exp(f(0))). Plugging then
 >> your new value as the old value into the right hand side gives you
 >> your next value, and so on.
 >> It may not be the most efficient way of integrating the equation, but
 >> it is probably the best way to keep control over your results as there
 >> is only one parameter required, namely dt. Just by making the latter
 >> smaller, it is easy to check if your solution has sufficiently
 >> converged or not.
 >>
 >> Thomas
 >>
 >
 >Thanks a lot Thomas. Are you suggesting a home-made approach to solve the 
 >ODE?
 >
 >If it works it would be great! My headache is that I got something out from 
 >Matlab's ODE45, and ODE 23, and they all give warning of dividing by zero 
 >and the results have NANs... I am not sure why... and want another approach 
 >which is guaranteed to be numerically stable and correct so I could compare 
 >the results ... 
 >
 > 

 only devel knows why you make the problem harder on every try. 
 your original problem was
  y'=c1*y*log(y)+c2*y+c3*y^r  
 we concluded (in agreement with Maple and Mathematica) that there is no 
 closed form solution. what you could do is separation of variables 

  integral{y0 to y} dy/(c1*y*log(y)+c2*y+c3*y^r) = integral{x0 to x} dx =x-x0
 and solve this pointwise for a given x in terms of the y. this involves a 
 quadrature routine and a nonlinear equation solver (for afficionatos of matlab:
 quadl and fzero) but this is quite costly. hence I would try a numerical 
 initil value solver, where you have to plug in nothing than (again in matlab)
 
 function yp=fun(x,y,param)  % param contains c1,c2,c3,r
 c1=param(1);
 c2=param(2);
 c3=param(3);
 r=param(4);
 if y <= 0 
  error(['integrator left domain of definition: y<=0 at x= ',num2str(x)]);
 end
 yp=c1*y*log(y)+c2*y+c3*y^r;

using this with ode45 for example you must also specify t0 and y0, the initial 
values. cleary this must not be t0=0,y0=0, otherwise you need special treatment 
of this since 0*log(0) will not evaluate to 0 but to NaN.
also, if you have some y0>0, it is not clear whether the solution exists 
for all x>x0, it might well end in a singularity and a crude grid might lead
to a discretized yi <0! 
hence you must take small steps, a decent integrator (ode15s (?)) 
and be aware of more numerical trouble (maybe also blowup).
hth
peter

On 22 May, 05:39, Robert Israel <isr...@math.MyUniversitysInitials.ca>
wrote:

> Without knowing the precise function f, I would suspect that either the
> solution is blowing up (y -> infinity) at a finite value of t, namely
> t = int_0^infty 1/(c1 y + c2 + c3 exp(f(y))) dy, or it's hitting a value
> of y where f is undefined.
>
> For example, try y' = y + 1 + exp(y), y(0) = 0; the solution goes to
> infinity at approximately t= .5871365611.
> --
> Robert Israel              isr...@math.MyUniversitysInitials.ca
> Department of Mathematics        http://www.math.ubc.ca/~israel
> University of British Columbia            Vancouver, BC, Canada

I get the same result for y' = y + 1 + exp(y) using my program that I
gave above, but I can't see that this should be due to a singularity
at a finite value of x. It is just so that the exp(y) term blows up
the solution so rapidly that it can't be represented anymore on the
computer. Mathematically it should still be finite (it only diverges
at infinity).

Thomas

Thomas Smid <thomas.smid@gmail.com> writes:

> On 22 May, 05:39, Robert Israel <isr...@math.MyUniversitysInitials.ca>
> wrote:
> 
> > Without knowing the precise function f, I would suspect that either the
> > solution is blowing up (y -> infinity) at a finite value of t, namely
> > t = int_0^infty 1/(c1 y + c2 + c3 exp(f(y))) dy, or it's hitting a value
> > of y where f is undefined.
> >
> > For example, try y' = y + 1 + exp(y), y(0) = 0; the solution goes to
> > infinity at approximately t= .5871365611.
> > --
> > Robert Israel              isr...@math.MyUniversitysInitials.ca
> > Department of Mathematics        http://www.math.ubc.ca/~israel
> > University of British Columbia            Vancouver, BC, Canada
> 
> I get the same result for y' = y + 1 + exp(y) using my program that I
> gave above, but I can't see that this should be due to a singularity
> at a finite value of x. It is just so that the exp(y) term blows up
> the solution so rapidly that it can't be represented anymore on the
> computer. Mathematically it should still be finite (it only diverges
> at infinity).

The solution is given implicitly by 

x = int_0^y(x) 1/(s + 1 + exp(s)) ds

As y -> infinity, x -> int_0^infty 1/(s + 1 + exp(s)) ds, which
is finite.  So the solution y(x) diverges as x approaches a finite
value.  

This situation is actually very common for nonlinear differential
equations.  If F(y) > 0 for y >= y_0 and 
lim inf_{y -> infty} ln(F(y))/ln(y) > 1, the solution of the ivp
y' = F(y), y(0) = y_0, will diverge as x approaches a finite value.
-- 
Robert Israel              israel@math.MyUniversitysInitials.ca
Department of Mathematics        http://www.math.ubc.ca/~israel 
University of British Columbia            Vancouver, BC, Canada

isr...@math.MyUniversitysInitials.ca> wrote:
> The solution is given implicitly by
>
> x = int_0^y(x) 1/(s + 1 + exp(s)) ds
>
> As y -> infinity, x -> int_0^infty 1/(s + 1 + exp(s)) ds, which
> is finite.  So the solution y(x) diverges as x approaches a finite
> value.
>
> This situation is actually very common for nonlinear differential
> equations.  If F(y) > 0 for y >= y_0 and
> lim inf_{y -> infty} ln(F(y))/ln(y) > 1, the solution of the ivp
> y' = F(y), y(0) = y_0, will diverge as x approaches a finite value.

Yes, I think you are right after all. Any power higher than 1 in y(t)
on the right hand side of y'(t)=.... should make the solution diverge
for finite t already. Actually, y'(t)=y(t)*ln(y(t)) should still be OK
(the solution should be y(t)=exp(exp(t)), which only diverges at
infinity). So, Mike's original equation   y'=c1*y*log(y)+c2*y+c3*y^r
should actually be better behaved in general (at least as long as
r<=1).

One should also be aware that the starting value y(t=0) can be all
important for how the system develops. For instance with the log(y)
term, it is crucial whether y is smaller or larger than 1 as this
determines the sign of log(y).

Anyway, if the differential equation has anything to do with a
realistic system, then the solution shouldn't diverge at all but stay
finite throughout the considered range. It is all a matter of
formulating a meaningful differential equation in the first place.
Then there shouldn't be any problems of this kind.

Thomas