Numerically integrating solution obtained from NDSolve method - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-11-18T08:19:49Z https://mathematica.stackexchange.com/feeds/question/58807 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/58807 1 Numerically integrating solution obtained from NDSolve method alekhine https://mathematica.stackexchange.com/users/8233 2014-09-03T19:17:46Z 2014-09-04T00:48:10Z <p>In the following example, \$u(x)\$ is found numerically using NDSolve method. </p> <pre><code> F = 1/1000 h = 12000/1000 d = 10/10 L = 1000 W = 3 phi[x_] := Piecewise[{{(1/2)*(1 - Tanh[((L*x)/(d))]), x &lt;= 1/2}, {(1/2)*(1 + Tanh[((L*(x - L/L))/(d))]), x &gt; 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 -&gt; "StiffnessSwitching", WorkingPrecision -&gt; 40, InterpolationOrder -&gt; 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 -&gt; All] </code></pre> <p>Which works perfectly. I need to see what is mean square error between obtained solution and another function \$vE(x)\$. </p> <pre><code>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] </code></pre> <p>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,</p> <pre><code>intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x, W*d/L, 1 - W*d/L}] </code></pre> <p>but this gives long error message ending with,</p> <pre><code> "...is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing" </code></pre> <p>How can I evaluate this integral?</p> https://mathematica.stackexchange.com/questions/58807/-/58830#58830 1 Answer by Michael E2 for Numerically integrating solution obtained from NDSolve method Michael E2 https://mathematica.stackexchange.com/users/4999 2014-09-04T00:48:10Z 2014-09-04T00:48:10Z <p>You can integrate the mean square error <code>mse</code> at the same time as computing <code>u[x]</code>.</p> <pre><code>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 -&gt; "StiffnessSwitching", WorkingPrecision -&gt; 40, InterpolationOrder -&gt; All]; Plot[Evaluate[mse[x] /. s], {x, W*d/L, 1 - W*d/L}, PlotRange -&gt; All] </code></pre> <p><img src="https://i.stack.imgur.com/TySWb.png" alt="Mathematica graphics"></p> <p>You seem to be interested in this change:</p> <pre><code>mse[1 - W*d/L] - mse[W*d/L] /. First@s (* 8198.070964448656291179158359833465311096 *) </code></pre> <hr> <p>The problem with the code</p> <pre><code> intVal = NIntegrate[({u[x]} /. s - vE[x])*({u[x]} /. s - vE[x]), {x, W*d/L, 1 - W*d/L}] </code></pre> <p>is a syntax issue. The expression that <code>/.</code> tries to apply is <code>s - vE[x]</code>. This can be seen from the <a href="http://reference.wolfram.com/mathematica/ref/TreeForm.html" rel="nofollow noreferrer"><code>TreeForm</code></a> of the expression:</p> <pre><code>Hold[({u[x]} /. s - vE[x])] // TreeForm </code></pre> <p><img src="https://i.stack.imgur.com/siJks.png" alt="Mathematica graphics"></p> <p>In other words, <code>ReplaceAll</code> has lower precedence than <code>Plus</code> (represented by the minus sign). The proper code is</p> <pre><code>intVal = NIntegrate[((u[x] /. First@s) - vE[x])^2, {x, W*d/L, 1 - W*d/L}, WorkingPrecision -&gt; 40] (* 8198.070964448656291179224357022321689725 *) </code></pre>