Re: Integro-differential analog of Love's equation via power series
- To: mathgroup at smc.vnet.net
- Subject: [mg79843] Re: [mg79823] Integro-differential analog of Love's equation via power series
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Mon, 6 Aug 2007 03:54:03 -0400 (EDT)
- References: <25633792.1186373747971.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
The code below is incomplete or otherwise doesn't work (at this machine,
anyway). f is undefined when Solve is executed, so f appears in "sol".
Then, when f[x_] is defined, we get infinite recursion.
Bobby
On Sun, 05 Aug 2007 04:00:01 -0500, chuck009 <dmilioto at comcast.com> wrote:
> 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];
>
>