Skip to main content

differential equations - 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?



Answer



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 *)

Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more