# Numerically integrating solution obtained from NDSolve method

In the following example, $u(x)$ is found numerically using NDSolve method.

 F = 1/1000
h = 12000/1000
d = 10/10
L = 1000
W = 3
phi[x_] :=
Piecewise[{{(1/2)*(1 - Tanh[((L*x)/(d))]),
x <= 1/2}, {(1/2)*(1 + Tanh[((L*(x - L/L))/(d))]), x > 1/2}}]
vE[x_] := x*(1 - x)*4
s = NDSolve[{u''[x] == (h*L*L/(d*d))*phi[x]*phi[x]*u[x] -
F*L*L*(1 - phi[x]), u[-W*d/L] == 0, u[1 + W*d/L] == 0},
u, {x, -W*d/L, 1 + W*d/L}, Method -> "StiffnessSwitching",
WorkingPrecision -> 40, InterpolationOrder -> All]
diff[x_] := (u[x] - vE[x])*(u[x] - vE[x])
Plot[Evaluate[{diff[x]} /. s], {x, W*d/L, 1 - W*d/L},
PlotRange -> All]


Which works perfectly. I need to see what is mean square error between obtained solution and another function $vE(x)$.

sum = 0;
Do[
first = W*d/L;
second = 1 - W*d/L;
{sum = sum + diff[first + (i/100)*(second - first)]},
{i, 0, 100, 1}]
Evaluate[sum]


but this gives only expression but not value. I think this is because $u(x)$ is obtained at discrete points only and is not defined on the points on which I have calculated error. I also tried using integration,

intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x,
W*d/L, 1 - W*d/L}]


but this gives long error message ending with,

 "...is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing"


How can I evaluate this integral?

• Evaluate[sum] gives number if you use diff[x_] := (u[x] - vE[x])*(u[x] - vE[x])/.s. Then NIntegrate[diff[x], {x, W*d/L, 1 - W*d/L}] works too. The current problem for your NIntegrate is that - has lower precedence than /., which you can easily check out by selecting the code click by click Sep 3, 2014 at 19:26
• on thing to note NDSolve returns a list of solutions, even though in this case there is only one solution, so you should be doing /. s[]. (as is you get a bunch of superfluous {} in the results ) Sep 3, 2014 at 20:45
• @Coolwater, it worked with your suggestion, thanks.. Sep 4, 2014 at 8:09

You can integrate the mean square error mse at the same time as computing u[x].

s = NDSolve[{
u''[x] == (h*L*L/(d*d))*phi[x]*phi[x]*u[x] - F*L*L*(1 - phi[x]),
u[-W*d/L] == 0, u[1 + W*d/L] == 0,
mse'[x] == (u[x] - vE[x])^2, mse[-W*d/L] == 0},
{u, mse}, {x, -W*d/L, 1 + W*d/L}, Method -> "StiffnessSwitching",
WorkingPrecision -> 40, InterpolationOrder -> All];

Plot[Evaluate[mse[x] /. s], {x, W*d/L, 1 - W*d/L}, PlotRange -> All] You seem to be interested in this change:

mse[1 - W*d/L] - mse[W*d/L] /. First@s

(* 8198.070964448656291179158359833465311096 *)


The problem with the code

 intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x, W*d/L, 1 - W*d/L}]


is a syntax issue. The expression that /. tries to apply is s - vE[x]. This can be seen from the TreeForm of the expression:

Hold[({u[x]} /. s - vE[x])] // TreeForm In other words, ReplaceAll has lower precedence than Plus (represented by the minus sign). The proper code is

intVal = NIntegrate[((u[x] /. First@s) - vE[x])^2, {x, W*d/L, 1 - W*d/L},
WorkingPrecision -> 40]

(* 8198.070964448656291179224357022321689725 *)

• Thanks a lot Michael for the complete explanation..!! Sep 4, 2014 at 8:08
• @alekhine You're welcome. :) Sep 4, 2014 at 11:41