1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Coolwater Sep 3 '14 at 19:26
  • $\begingroup$ 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[[1]]. (as is you get a bunch of superfluous {} in the results ) $\endgroup$ – george2079 Sep 3 '14 at 20:45
  • $\begingroup$ @Coolwater, it worked with your suggestion, thanks.. $\endgroup$ – alekhine Sep 4 '14 at 8:09
1
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

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 *)
$\endgroup$
  • $\begingroup$ Thanks a lot Michael for the complete explanation..!! $\endgroup$ – alekhine Sep 4 '14 at 8:08
  • $\begingroup$ @alekhine You're welcome. :) $\endgroup$ – Michael E2 Sep 4 '14 at 11:41

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.