Solving an algebraic-integral equation in Mathematica

[nickthequick]

Hi,

I'm trying to solve an algebraic-integral equation (I don't know if this is the proper description of this class of equations, it just seems like the least wrong way to describe them) and have run into several issues that I'll describe below, but first I'll outline the problem.

I'm trying to solve for [itex]k[/itex] dependent variables where the governing equation for the [itex]k-th[/itex] mode takes the form

[tex] \alpha_{k} +\frac{1}{\lambda}\gamma \int_0^{\lambda} \frac{(1+\mathcal{X})\cos k\xi +\mathcal{Y}\sin k\xi }{\sqrt{(1+\mathcal{X})^2+\mathcal{Y}^2} } d\xi =0 [/tex]

where

[tex] \mathcal{X}=2\sum_{n=1}^N n(\alpha_n \cos n\xi);\quad \& \quad \quad \mathcal{Y} = 2\sum_{n=1}^N n(-\alpha_n \sin n\xi ).
[/tex]

and the parameters [itex]\lambda, \gamma [/itex] and N are prescribed, where N is a parameter that will be be increased to meet convergence criterion. The independent variable of the system is [itex]\gamma[/itex] and we would like to solve for the [itex]\alpha_k[/itex].

Now, I don't really know how to deal with the integral term and have no reason to suspect there's a better way to rewrite it, so I am doing the most naive thing that comes to mind, namely, taking the so called trapezoidal approximation of the integral over even spacing. That is, given a function [itex]f[/itex] on an interval [itex](a,b)[/itex], we have

[tex]\int_a^b f(x) dx =\frac{b-a}{2M} \left(f(x_1)+2\sum_{i=2}^Mf(x_i) +f(x_{M+1})\right) +R_M[/tex]

where [itex]x_i=a+\frac{b-a}{M}*i [/itex].

My goal is to increase the resolution of the integral until the remainder [itex] R_M[/itex] is sufficiently small.

I'm having a few problems. My main question is, is this the proper way to go about problems of this nature? Something is wrong with either my physics (which I understand you cannot help with without more details, which I can provide if interested) or my code as I get very strange non-convergent behavior as M is increase. Secondly, I am finding that for M<N, the gross majority of the energy tends to congregate in the M-th mode. I'm not completely sure why the system would bias this mode so I've just been considering cases where M>N. Last, but certainly not least, this code ends up being very memory intensive for test examples that I'm playing with, so I'm looking at ways of getting around this.

Any suggestions for my code, or even on how to approach these types of problems, would be much appreciated.


My mathematica code is included below.
M = 20;
\[Kappa] = 10^6;
\[Lambda] = 2 \[Pi]/\[Kappa];
\[Gamma] = 1.0126;
f1a[k_, x_] := ((1 +
2 Sum[n*(\[Alpha][n]*Cos[2 \[Pi]*n*x]), {n, 1, M}]) Cos[
2 \[Pi]*k*x] +
2*Sin[2 \[Pi]*k*x]*
Sum[n*(-\[Alpha][n]*Sin[2 \[Pi]*n*x]), {n, 1, M}])/
Sqrt[(1 + 2 Sum[n*(\[Alpha][n]*Cos[2 \[Pi]*n*x]), {n, 1, M}])^2 +
Sum[2 n*(-\[Alpha][n]*Sin[2 \[Pi]*n*x]), {n, 1, M}]^2];
MMax = 22;
I1a[j_] := 1/(MMax) (Sum[f1a[j, (i - 1)/(MMax)], {i, 1, MMax}]);
Governing[j_] := \[Alpha][j] - \[Gamma]*I1a[j];
IC = Join[Table[{\[Alpha][1], .4}, {i, 1}],
Table[{\[Alpha], 0}, {i, 2, M}]];
S =
FindRoot[Table[With[{i = i}, Governing == 0], {i, 1, M}],IC];
Y[x_] := Evaluate[
Sum[-2*(\[Alpha][k]^2)/\[Kappa] +
2*Cos[2 \[Pi]*x*k]*\[Alpha][k]/\[Kappa], {k, 1, M}] /.
S];
X2[x_] :=
x*\[Lambda] +
Evaluate[
Sum[Sin[2 \[Pi]*x*k]*\[Alpha][k]/\[Kappa], {k, 1, M}] /. S];
ParametricPlot[{X2[x], Y[x]}, {x, 0, 1}, AspectRatio -> 1/4]


Nick

 

[SteamKing]

Integral equations can be tricky to solve, and the numerical methods which can solve them involve using a lot of complicated mathematics to figure out the solutions.

Because of the nature of the integral expression, the use of the trapezoidal method to calculate the integral may not provide accurate results. There are numerical integration methods which are specially adapted to give accurate results when trigonometric or other periodic functions are involved.