MathGroup Archive 2007

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

Search the Archive

Integro-differential analog of Love's equation via power series

  • To: mathgroup at smc.vnet.net
  • Subject: [mg79823] Integro-differential analog of Love's equation via power series
  • From: chuck009 <dmilioto at comcast.com>
  • Date: Sun, 5 Aug 2007 05:00:01 -0400 (EDT)

Hello guys,

I modified Paul Abbott's code (you mind Paul?) to solve the integro-differential analog of Love's equation:

y'[x]=1+1/pi Integrate[y[t]/((x-t)^2+1),{t,-1,1}]

using a simple power series:

1.  I have to specify an initial condition y0=1/2=c0

2.  The left side now contains the derivative of the power series.

3.  The right side is similar to the integral equation except that now I'm solving for just 
c1 through cn although I keep c0 in the matrix for now unless someone can suggest a better way.

4.  When I use Solve, I replace c0 with 1/2 and also specify that now I'm only solving for {c1, c2,...cn}.  I use Drop to remove c0 from the list.

5.  For the check I explicityly calculate the derivative of the resulting power series and then plot the difference lhs-rhs of the IDE.  The result with 40 equations and 40 unknowns is less than 10^{-6} execpt for spikes at the end-points and also receive a "row-reduce of badly-conditioned matrix".

So, anyone wants to look at it and suggest more concise ways of coding it or explain why I would obtain a badly-conditioned matrix?

a = -1; 
b = 1; 
y0 = 1/2; 
n = 40;

xs = N[Table[x, {x, -1, 1, (b - a)/(n - 1)}], 6]; 
cs = Table[Subscript[c, k], {k, 0, n}]; 

lhs = cs . Table[(i - 1)*xs^(i - 2), {i, 1, n + 1}]; 

rhs = f[xs] + (1/Pi)*cs . Table[NIntegrate[
        Evaluate[t^i/((xs - t)^2 + 1)], {t, -1, 1}], {i, 0, n}]; 

sol = Solve[lhs == rhs /. Subscript[c, 0] -> y0, Drop[cs, 1]]

f[x_] = Sum[Subscript[c, i]*x^i, {i, 0, n}] /. First[sol] /. 
   Subscript[c, 0] -> y0

Plot[f[x], {x, -1, 1}]
fd[x_] = D[f[x], x]
Plot[fd[x], {x, -1, 1}]

Plot[1 + (1/Pi)*NIntegrate[f[t]/((x - t)^2 + 1), {t, -1, 1}] - 
    fd[x], {x, -1, 1}, PlotRange -> All];


  • Prev by Date: Re: Documentation Center (v6): do-it-yourself Mathematica
  • Next by Date: Re: ListPlot3D on 3 by 3 Arrays
  • Previous by thread: Re: Re: request for a few minutes CPU-time
  • Next by thread: Re: Integro-differential analog of Love's equation via power series