1
$\begingroup$

We consider the following integral equation of mixed type: $x[t]=\frac{t^4}{6} -\frac{t^3}{3} +t+\Sigma_{\mu=1}^2\int_{0}^t k_\mu[t,s]G_\mu[s,x[s]]ds ,0\leq t\leq 1 $ ; (1)
where

$k_1 [t, s] = s^3$

$k_2 [t, s] = -2 (t - s)$

$G_1 [t,x[t]]=\frac{1}{x[t]}$

$G_2 [t,x[t]]=x^2 [t]$

1.Let us the basis functions are piecewise linear function and the collocation points are the zeros of Chebychev polynomials. Our aim is to solve the integral equation (1) by the standard collocation method for $N=6,8$
In this part, we have problem to solve integral because the denominator is zero. We may need to approximate integrand by taylor expand finite, but how to do that?

Clear["Global`*"]
T[1, t_] = t;
T[2, t_] = (4 t^2) - 1;
T[n_, t_] := 2 *t*T[n - 1, t] - T[n - 2, t];
tableoft = tt /. NSolve[T[6, tt], tt];
Print["The roots of 3th Chebyshev polynomial (or collocation points) \
is:", "\n", tableoft];

ttt = (1/2) (tableoft + 1)
t[j_] := j/6
u[i_, kt_] := Which[
  i == 1, If[t[1] <= kt <= t[2], (kt - t[2])/(t[1] - t[2]), 0],
  i == 6, u[1, 1 - kt],
  i >= 2 && i < 6, 
  Switch[kt, t[i - 1] < kt <= t[i], (kt - t[i - 1])/(t[i] - t[i - 1]),
    t[i] < kt <= t[i + 1], (kt - t[i + 1])/(t[i] - t[i + 1]); 0]]
x[t_] := Sum[Simplify[Subscript[a, i]*u[i, t]], {i, 1, 6}]
f[t_] := (t^4/6 - t^3/3 + t)
Subscript[k, 1][t_, s_] := s^3;
Subscript[G, 1][t_] := 1/x[t];
Subscript[k, 2][t_, s_] := -2 (t - s)
Subscript[G, 2][t_] := (x[t])^2;
pp[t_, s_] := Subscript[k, 1][t, s]*Subscript[G, 1][s] // N;
ppp[t_, s_] := Subscript[k, 2][t, s]*Subscript[G, 2][s] // N;
pp1[t_] := Simplify[Integrate[pp[t, s], {s, 0, t}]];
pp2[t_] := Simplify[Integrate[ppp[t, s], {s, 0, t}]];
d = Table[
   Expand[(x[ttt[l]] - f[ttt[[l]]] + pp1[ttt[[l]]] + 
       pp2[ttt[[l]]])] == 0, {l, 1, 6}];
list = Flatten[Table[{Subscript[a, i], Subscript[b, i]}, {i, 1, 6}]]
NSolve[d, list]

2.Now let $z_r [t]=G_r [t,x[t]] , r=1,2 $ thus the form of the function $z_r [t]$ is:
$z_r [t]=G_r [t,\frac{t^4}{6} -\frac{t^3}{3} +t+\Sigma_{\mu=1}^2\int_{0}^t k_\mu[t,s]z_\mu[s]ds , r=1,2 ; (2)$
i.e, we obtain a system of nonlinear integral equations. Suppose that its solution is:$z[t]=(z_1 [t],z_2 [t])$,then the solution of (1) is given by
$x[t]=\frac{t^4}{6} -\frac{t^3}{3} +t+\Sigma_{\mu=1}^2\int_{0}^t k_\mu[t,s]z_\mu[s]ds$
where $z_r [t]$ is the solution of the system of integral equation (2).
I try to solve the implicit nonlinear mixed integral equation (2) using the collocation method with the same conditions as the first case.
Notice that piecewise linear functions $u_i [t]$ are defined in code.

I don't know why we got to false equation in ddd[[1]]!

I attached below the system of 12 equations that I am not able to solve it

Clear["Global`*"]

T[1, t_] = t;

T[2, t_] = (4 t^2) - 1;

T[n_, t_] := 2 *t*T[n - 1, t] - T[n - 2, t];

tableoft = tt /. NSolve[T[6, tt], tt];

Print["The roots of 6th Chebyshev polynomial (or collocation points) \

is:", "\n", tableoft];

ttt = (1/2) (tableoft + 1)

t[j_] := j/6
u[i_, kt_] := Which[
   i == 1, 
   Piecewise[{{(kt - t[2])/(t[1] - t[2]), t[1] <= kt <= t[2]}}, 0],
   i == 6, u[1, 1 - kt],
   i >= 2 && i < 6, 
   Piecewise[{{(kt - t[(i - 1)])/(t[i] - t[(i - 1)]), 
      t[i - 1] < kt <= t[i]}, {(kt - t[(i + 1)])/(t[i] - t[(i + 1)]), 
      t[i] < kt <= t[(i + 1)]}}, 0]];

f[t_] := (t^4/6 - t^3/3 + t);

Subscript[k, 1][t_, s_] := s^3;

Subscript[k, 2][t_, s_] := -2 (t - s);

S1[t_] := 
  Simplify[Integrate[
    Subscript[k, 1][t, s]*Sum[Subscript[a, i]*u[i, s], {i, 1, 6}], {s,
      0, t}]];

S2[t_] := 
  Simplify[Integrate[
    Subscript[k, 2][t, s]*Sum[Subscript[b, i]*u[i, s], {i, 1, 6}], {s,
      0, t}]];
d = Table[
   Expand[(S1[ttt[[l]]] + S2[ttt[[l]]] + f[ttt[[l]]])*
       Sum[Subscript[a, i]*u[i, ttt[[l]]], {i, 1, 6}] - 1 == 0], {l, 
    1, 6}];
dd = Table[
   Expand[Sum[
       Subscript[b, i]*u[i, ttt[[l]]], {i, 1, 
        6}] - ((S1[ttt[[l]]] + S2[ttt[[l]]] + f[ttt[[l]]])^2) == 
     0], {l, 1, 6}];
    ddd = Flatten[{d, dd}];
list = Flatten[Table[{Subscript[a, i], Subscript[b, i]}, {i, 1, 6}]];

NSolve[ddd, list]

Mathematica returned below answer

{}

Notice that 2 equations obtained false! Is this a problem for solve this 12 * 12 equations?

$\endgroup$
4
  • $\begingroup$ ddd[[1]] is False. $\endgroup$
    – cvgmt
    Jul 23 at 9:10
  • 1
    $\begingroup$ If any equations evaluate to False, you can never solve the system. There is no value of x that solves 1 == 2. $\endgroup$ Jul 23 at 10:13
  • 1
    $\begingroup$ In the first collocation point Sum[Subscript[a, i]*u[i, ttt[[1]]], {i, 1, 6}]=0, it is why d[[1]] is False. $\endgroup$ Jul 23 at 10:40
  • $\begingroup$ Technically you try to solve integral equation with using wavelets based collocation method, and not 12 nonlinear equations. $\endgroup$ Jul 23 at 19:22
3
$\begingroup$

I can recommend collocation method based on Bernoulli wavelets described here. With this method we can solve problem for 8 and 16 collocation points (n=2 or n=3 in this code):

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"]; \
Get["NumericalDifferentialEquationAnalysis`"];

ue[t_] := t;
f[t_] := t - t^3/3 + t^4/6;
n = 3;
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1;
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}]; xcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}];
psi1[x_] := Piecewise[{{BernoulliB[2, x], 0 <= x < 1}, {0, True}}];
psi2[x_] := Piecewise[{{BernoulliB[1, x], 0 <= x < 1}, {0, True}}];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk[x_, j_, k_] := (psi1jk[x, j, k] + psi2jk[x, j, k])/2;


np = 2 M; g = GaussianQuadratureWeights[np, -1, 1, 60]; points = 
 g[[All, 1]];
weights = g[[All, 2]];
Int[ff_, z_] := Sum[(ff /. z -> points[[i]])*weights[[i]], {i, 1, np}];
u[t_] := Sum[
    a[j, k]*psijk[t, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}] + a0;
int[t_] := 
  t/2 Int[(t/2 (1 + z))^3/u[t/2 (1 + z)] - 
     2 (t - t/2 (1 + z)) u[t/2 (1 + z)]^2, z](*s\[Rule]x/2 (1+z)*);
eq = Table[-u[xcol[[i]]] + f[xcol[[i]]] + int[xcol[[i]]] == 0, {i, 
    Length[xcol]}];
varM = Join[{a0}, 
   Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];

sol = FindRoot[eq, Table[{varM[[i]], 1/10}, {i, Length[varM]}], 
   WorkingPrecision -> 30, MaxIterations -> 1000];
unum = Table[{xcol[[i]], Evaluate[u[xcol[[i]]] /. sol]}, {i, 
    Length[xcol]}];

du = Table[{x, Abs[ue[x] - Evaluate[u[x] /. sol]]}, {x, xcol}]; 

Visualization for 8 (upper line) and 16 collocation points

{Show[Plot[ue[x], {x, 0, 1}, 
   PlotLegends -> 
    Placed[LineLegend[{"Exact"}, LabelStyle -> {Black, 15}], 
     Scaled[{0.7, 0.2}]], AspectRatio -> 1, 
   LabelStyle -> Directive[{FontSize -> 15}, Black], 
   AxesLabel -> {"t", "x"}, PlotStyle -> Blue], 
  ListPlot[unum, PlotRange -> All, PlotStyle -> Red, 
   PlotLegends -> 
    Placed[PointLegend[{"Numeric"}, LabelStyle -> {Black, 15}], 
     Scaled[{0.7, 0.3}]]]], 
 ListPlot[du, Filling -> Axis, PlotLabel -> "Absolute error", 
  PlotRange -> All]}

Figure 1

With some modification we can also solve integral equation with 6 colocation points and with absolute error of $5.5\times 10^{-17}$ as follows

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"]; \
Get["NumericalDifferentialEquationAnalysis`"]; Clear["Global`*"]

k0 = 2; M0 = 3;

T[m_, t_] := BernoulliB[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) Sqrt[2/Pi] T[m, 2^k t - 2 n + 1], (n - 1)/
       2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
Psi[k_, M_, t_] := 
  Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]];
With[{k = k0, M = M0}, 
 var = Flatten[Table[a[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]; 
x[t_] := var . Psi[k0, M0, t];(*tableoft=tt/. NSolve[T[8,tt],tt];
ttt=(1/2) (tableoft+1);*)

dt = 1/Length[var]; tl = Table[ l*dt, {l, 0, Length[var]}]; ttt = 
 Table[(tl[[l - 1]] + tl[[l]])/2, {l, 2, Length[var] + 1}]; np = 
 2 Length[var]; g = GaussianQuadratureWeights[np, -1, 1, 60]; points =
  g[[All, 1]];
weights = g[[All, 2]];
Int[ff_, z_] := Sum[(ff /. z -> points[[i]])*weights[[i]], {i, 1, np}];


f[t_] := t^4/6 - t^3/3 + t;
k1[t_, s_] := s^3;
G1[s_] := 1/x[s];
k2[t_, s_] := -2 (t - s)
G2[s_] := x[s]^2;
pp[t_, s_] := k1[t, s]*G1[s];
ppp[t_, s_] := k2[t, s]*G2[s];
pp1[t_] := 
 t/2 Sum[(pp[t, t/2 (1 + z)] /. z -> points[[i]])*weights[[i]], {i, 1,
     np}]; s = t/2 (1 + z);
pp2[t_] := 
 t/2 Sum[(ppp[t, t/2 (1 + z)] /. z -> points[[i]])*weights[[i]], {i, 
    1, np}]; 
int[t_] := 
 t/2 Int[(t/2 (1 + z))^3/x[t/2 (1 + z)] - 
    2 (t - t/2 (1 + z)) x[t/2 (1 + z)]^2, z]
d = Table[-x[t] + f[t] + int[t], {t, ttt}];
sol = FindRoot[Table[d[[i]] == 0, {i, 1, Length[d]}], 
   Table[{var[[i]], 1/10}, {i, Length[var]}], WorkingPrecision -> 30, 
   MaxIterations -> 200];

Visualization

Show[Plot[t, {t, 0, 1}, AxesLabel -> {"t", "x"}], 
 ListPlot[Table[ {ttt[[i]], Evaluate[x[ttt[[i]]] /. sol]}, {i, 
    Length[ttt]}], PlotStyle -> Red]]

Plot[Abs[t - x[t] /. sol], {t, 0, 1}, PlotRange -> All, 
 PlotLabel -> "Absolute Error"]

Figure 2

Finally we can repair last code @user68119. There are several typos with u[i,kt] definition, and after all we have

 Clear["Global`*"]
n = 8;
T[1, t_] = t;

T[2, t_] = (4 t^2) - 1;

T[n_, t_] := 2*t*T[n - 1, t] - T[n - 2, t];

tableoft = tt /. NSolve[T[n, tt], tt];

Print["The roots of 6th Chebyshev polynomial (or collocation points) 
is:", "\n", tableoft];

ttt = (1/2) (tableoft + 1)

t[j_] := ttt[[j]]; 
u[1, kt_] := 
 Piecewise[{{(kt - t[2])/(t[1] - t[2]), t[1] <= kt <= t[2]}}, 0]; 
u[n, kt_] := u[1, 1 - kt];
Do[u[i, kt_] := 
   Piecewise[{{(kt - t[(i - 1)])/(t[i] - t[(i - 1)]), 
      t[i - 1] < kt <= t[i]}, {(kt - t[(i + 1)])/(t[i] - t[(i + 1)]), 
      t[i] < kt <= t[(i + 1)]}}, 0], {i, 2, n - 1}];

f[t_] := (t^4/6 - t^3/3 + t);

Subscript[k, 1][t_, s_] := s^3;

Subscript[k, 2][t_, s_] := -2 (t - s);

S1[t_] := 
  Simplify[Integrate[
    Subscript[k, 1][t, s]*Sum[Subscript[a, i]*u[i, s], {i, 1, n}], {s,
      0, t}]];

S2[t_] := 
  Simplify[Integrate[
    Subscript[k, 2][t, s]*Sum[Subscript[b, i]*u[i, s], {i, 1, n}], {s,
      0, t}]];

d = Table[
   Expand[(S1[ttt[[l]]] + S2[ttt[[l]]] + f[ttt[[l]]])*
       Sum[Subscript[a, i]*u[i, ttt[[l]]], {i, 1, n}] - 1 == 0], {l, 
    1, n}];

dd = Table[
   Expand[Sum[
       Subscript[b, i]*u[i, ttt[[l]]], {i, 1, 
        n}] - ((S1[ttt[[l]]] + S2[ttt[[l]]] + f[ttt[[l]]])^2) == 
     0], {l, 1, n}];

ddd = Flatten[{d, dd}];
list = Flatten[Table[{Subscript[a, i], Subscript[b, i]}, {i, 1, n}]];

sol = FindRoot[ddd, Table[{list[[i]], 1/10}, {i, Length[list]}]]

Visualization for n=6 and n=8

    i1 = Evaluate[
  Subscript[k, 1][t, s]*
   Sum[Subscript[a, i]*u[i, s] /. sol, {i, 1, n}]]; i2 = 
 Evaluate[Subscript[k, 2][t, s]*
   Sum[Subscript[b, i]*u[i, s] /. sol, {i, 1, n}]];

x[t_?NumericQ] := 
 f[t] + Quiet[NIntegrate[i1, {s, 0, t}] + NIntegrate[i2, {s, 0, t}]]

xnum = Table[{t, x[t]}, {t, ttt}]

{Show[Plot[t, {t, 0, 1}, AxesLabel -> {"t", "x"}], 
  ListPlot[xnum, PlotStyle -> Red]], 
 ListPlot[Abs[ttt - xnum[[All, 2]]], Filling -> Axis, 
  PlotLabel -> "Absolute error"]}

Figure 3 Figure 4

$\endgroup$
1
$\begingroup$

The first element of d is equal to exactly -1. This is because all of the elements u[i, ttt[[1]]] are zero:

With[{l = 1}, Table[u[i, ttt[[l]]], {i, 1, 6}]]
(* {0,0,0,0,0,0} *)

For the same reason, the first element of dd is equal to -(S1[ttt[[1]]] + S2[ttt[[1]]] + f[ttt[[1]]])^2. It happens that S1[ttt[[1]]] and S2[ttt[[1]]] are zero as well, so this works out to f[ttt[[1]]]^2:

With[{l = 1}, {S1[ttt[[l]]], S2[ttt[[l]]], f[ttt[[l]]]}]
(* {0., 0., 0.0749211} *)

As neither -1 nor -f[ttt[[l]]]^2 are equal to zero, this means that the equations in which these are equated to zero automatically evaluate to False.

As to what these quantities should be equal to, I am not familiar enough with the mathematical method you describe to know for certain. But given the above, I would suspect that the definition of u is where your error is.

$\endgroup$
1
$\begingroup$

Not an answer for your desired collocation method, but extented comment.

I don't know whether you already know the solution of this integral equation. It is x[t] = t

Get this with two fold differentiation and then NDSolve.

eq[t_] = x[t] == 
t - t^3/3 + t^4/6 + Integrate[s^3 1/x[s] , {s, 0, t}] + 
Integrate[-2 (-s + t) x[s]^2 , {s, 0, t}]

eq[0]
(*   x[0] == 0   *)

D[eq[t], t] /. t -> 0
(*   Derivative[1][x][0] == 1   *)

d2eq = D[eq[t], {t, 2}]

xsol[a_, b_] := 
  x /. First@NDSolve[{d2eq, x[0] == a, x'[0] == b}, x, {t, 0, 1}]

Since a == 0 gives singularity

Plot[Evaluate[{xsol[10^-10, 1][t]}], {t, 0, 1}]

This implies x[t] = xsol[0,1] = t. Indeed

eq[t] /. x -> Function[t, t]

(*   True   *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.