1

I am trying to solve the following integral equation numerically in mathematica:

$f(x) = 1 + \int_0^1 dy \left(1-3xy\right)f(y)$. The exact solution is $f(x)=\frac{8}{3}-2x$. I can get the desired result for fix number of iteration as follows:

dim = 11;
<< NumericalDifferentialEquationAnalysis`;
dat = GaussianQuadratureWeights[dim, 0, 1];
dat2 = Table[dat[[i, 2]], {i, 1, dim}];
approxsoln = Table[1, {i, 1, dim}];
f = Table[1, {i, 1, dim}];
Kn = Table[(1 - 3 dat[[i, 1]] dat[[j, 1]]), {i, 1, dim}, {j, 1, 
dim}];
exactsoln = Table[(8/3 - 2 dat[[i, 1]]), {i, 1, dim}];

Do[values = f + Total[Kn approxsoln dat2, 1];
 Print[iter, ",  " , approxsoln[[1]] , ",  ", exactsoln[[1]]];
 approxsoln = values, {iter, 1, 20}]

Note that I am using Gaussian Quadrature for the numerical integration for the faster evaluation of the actual problem. Now, instead of fix number iteration, I am trying to set a fix accuracy to stop the iteration ($\chi^2$ test ) as follows: 

chsq=0.;
acc=0.001;
Do[values = f + Total[Kn approxsoln dat2, 1];
chsq = Max[Abs[1 - values/approxsoln], chsq]; 
If[chsq < acc, Break[]]; approxsoln = values; 
Print[iter, ",  " , approxsoln[[1]] , ",  ", exactsoln[[1]]], {iter, 
1, 20}]

But unable to stop the iteration. Does anybody have any idea what am I doing wrong or how to do this? Thanks in advance for the help and for your time.

H.N.
  • 21
  • 6
  • If all you want is to use the Gauss-Legendre abscissas and weights, there's no need to load that package. See this, for instance. – J. M.'s missing motivation Feb 01 '17 at 14:16
  • Thanks. Yes, I know I can define it locally instead of loading the package. As it is available, I am just loading that for instance. Additionally, that won't solve my problem. All I want is to stop the iteration at some fix accuracy. – H.N. Feb 01 '17 at 14:22
  • Probably better to use While or NestWhile than Do with an If-Break. – march Feb 01 '17 at 16:45
  • While playing with this: when I choose an initial value for chsq, it does not change with each iteration. You probably want to check that part of the code. – march Feb 01 '17 at 16:47
  • you should reexamine the logic of this expression: chsq = Max[Abs[1 - values/approxsoln], chsq]; If your solution improves you keep the last solution error measure (why?). – george2079 Feb 01 '17 at 18:46
  • 1
    the curly brackets here acc=10^{-3} are a problem too. This makes acc a list and so your logical test chsq < acc is never true as you compare a scalar to a list. – george2079 Feb 01 '17 at 18:50
  • Sorry about acc=10^{-3}. I just mistakenly wrote here in that form. – H.N. Feb 01 '17 at 18:59

1 Answers1

1

I want to thank all for your time and for you valuable suggestions. I find the answer by myself by changing the code a bit. I am posting it here.

dim = 11;
<< NumericalDifferentialEquationAnalysis`;
dat = GaussianQuadratureWeights[dim, 0, 1];
dat2 = Table[dat[[i, 2]], {i, 1, dim}];
approxsoln = Table[1, {i, 1, dim}];
f = Table[1, {i, 1, dim}];
Kn = Table[(1 - 3 dat[[i, 1]] dat[[j, 1]]), {i, 1, dim}, {j, 1, dim}];
exactsoln = Table[(8/3 - 2 dat[[i, 1]]), {i, 1, dim}];

Do[
values = f + Total[Kn approxsoln dat2, 1];
chsq = Max[Abs[1 - values/approxsoln]];
If[chsq < .000001, Break[]];
Print[iter, ",  " , approxsoln[[1]] , ",  ", exactsoln[[1]]];
approxsoln = values,{iter,1,60}]
H.N.
  • 21
  • 6