Skip to main content

differential equations - NDSolve in NIntegrate function, numerical values as a symbolic variables


I'm trying to implement the following iterative scheme $g_{n+1}=g_n-\int_0^t\Phi(g_n) dt$, where $g_n$ is a source term in a PDE and $\Phi$ is the solution of an adjoit problem associated to the PDE which use the solution of the PDE as final data. to this aim we need many steps :



0) Take $g_0=0$.


1) Solve the PDE to obtain the solution $u(t,x,g_n)$.


2) Solve the adjoint PDE to obtain $\Phi(t,x,g_n)$.


3) Calculate $g_{n+1}(x)=g_n(x)-\int_0^t\Phi(t,x,g_n) dt$, and go to 1).




I tried this one


nsol = NDSolve[{D[u[t, x], t] == D[u[t, x], x, x] + #, u[0, x] == 1, 
u[t, 0] == u[t, 1] == 1}, u, {t, 0, 1}, {x, 0, 1}] &; (* Solve the pde with source # *)

nasol = NDSolve[{D[v[t, x], t] == -D[v[t, x], x, x], 
 v[1, x] == First[u[1, x] /. nsol[#]] - umes[x], 
 v[t, 0] == v[t, 1] == -0.05}, v, {t, 0, 1}, {x, 0, 1}] &; (* Solve the adjont problem *)

Phi := NIntegrate[Evaluate[First[v[t, #2] /. nasol[#1[x]]]], {t, 0, 1}]&;

(* Calculate the integral in step 3). Here is the problem !!*)

g[0][x_] := 0;
umes[x_] := First[u[1, x] /. nsol[g[0][x]]] + 0.05;
g[n_Integer?Positive][x_] := g[n - 1][x] - Phi[g[n - 1], x];(* iteration *)

Now when I test the program, it seems to work in a good way


In[68]:= Phi[g[0], 0.1]
Out[68]= -0.05
In[70]:= g[1][0.1]

Out[70]= 0.05

To turn the loop for calculating the function $g_2$ we need to put $g_1$ in step 1) and here is the problem : The functions Phi[g[1],x] is just a numerical function and not working for symbolic variable $x$ because of the NIntegrate. Hence, the same problem for g[1][x].


Here is the Mathematica error


In[80]:= Phi[g[0], x]


NIntegrate::inumr: The integrand <<1>> has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.0547861}}.



I tried to use NumericQ for the variables but still not working. Also, I don't get how to close the loop using stopping condition like Norm[Phi(g[n][x],x]]<10^{-6}.



Finally, I notice that I'm new to Mathematica. Thanks for any help.



Answer



I show an algorithm for solving such problems. Calculation back is eliminated by replacing t'->-t, {-1,0}->{0,1}:


g[0][x_] := .01*x; n = 5;

Do[nsol[i] = 
  NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x] + g[i - 1][x], 
    u[0, x] == 1, u[t, 0] == 1, u[t, 1] == 1}, 
   u, {t, 0, 1}, {x, 0, 1}];
 nasol[i] = 

  NDSolveValue[{D[v[t, x], t] == D[v[t, x], x, x], 
    v[0, x] == nsol[i][1, x] - nsol[1][1, x] - .05, v[t, 0] == -.05, 
    v[t, 1] == -0.05}, v, {t, 0, 1}, {x, 0, 1}];
 g[i] = Interpolation[
   Table[{x, g[i - 1][x] - NIntegrate[nasol[i][t, x], {t, 0, 1}]}, {x,
      0, 1, .1}]];, {i, 1, n}]



 {Plot[Evaluate[Table[g[i][x], {i, 0, n}]], {x, 0, 1}, 

  AxesLabel -> {"x", "g"}], 
 Plot3D[nsol[n][t, x], {t, 0, 1}, {x, 0, 1}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> {"t", "x", ""}, 
  PlotLabel -> "nsol[n]"], 
 Plot3D[nasol[n][t, x], {t, 0, 1}, {x, 0, 1}, Mesh -> None, 
  ColorFunction -> Hue, AxesLabel -> {"t", "x", ""}, 
  PlotLabel -> "nasol[n]", PlotRange -> All]}

fig1


In the case of 2D +1, we use summation instead of the NIntegrate[]. In this example npoints=8.



<< NumericalDifferentialEquationAnalysis`
gl[npoints_] := 
 Block[{npo = npoints}, {pts, w} = 
   Transpose[GaussianQuadratureWeights[npo, 0, 1]]; {w, pts, npo}]

f[x_, y_] := 1; (* The exact source term to be constructed *)
Plot3D[f[t, x], {t, 0, 1}, {x, 0, 1}, 
 ColorFunction -> "TemperatureMap", AxesLabel -> {"t", "x", ""}, 
 PlotLabel -> "Source term f(x,y)", PlotLegends -> Automatic]


nsoleq = NDSolveValue[{D[u[t, x, y], t] == 
     D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + f[x, y], 
    u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0, 
    u[t, x, 1] == 0, u[t, 1, y] == 0}, 
   u, {t, 0, 1}, {x, 0, 1}, {y, 0, 
    1}]; (* nsoleq[1,x,y] is the observation used in construction of \
the source term *)
Plot3D[nsoleq[1, x, y], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "TemperatureMap", AxesLabel -> {"x", "y", ""}, 
 PlotLabel -> "t=1", PlotLegends -> Automatic]

g[0][x_, y_] := 0.;  (* Initialization of the iteration *)



With[{np0 = 16, np1 = .05, n = 40}, 
 Do[nsol[i] = 
   NDSolveValue[{D[u[t, x, y], t] == 
      D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + g[i - 1][x, y], 
     u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0, 
     u[t, x, 1] == 0, u[t, 1, y] == 0}, 

    u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}];

  nasol[i] = 
   NDSolveValue[{D[v[t, x, y], t] == 
      D[v[t, x, y], x, x] + D[v[t, x, y], y, y], 
     v[0, x, y] == nsol[i][1, x, y] - nsoleq[1, x, y], 
     v[t, x, 0] == 0, v[t, 0, y] == 0, v[t, x, 1] == 0, 
     v[t, 1, y] == 0}, v, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 

        "MinPoints" -> 5*15 + 1, "MaxPoints" -> 5*15 + 1, 
        "DifferenceOrder" -> Automatic}}];
  pp = Interpolation[
    Chop[Flatten[
      Table[{{x, y}, 
        Sum[nasol[i][gl[np0][[2]][[j]], x, y]*gl[np0][[1]][[j]], {j, 
          1, np0}]}, {x, 0, 1, np1}, {y, 0, 1, np1}], 1]]];
  p[i][x_, y_] := pp[x, y] + 0.000001*g[i - 1][x, y];
  nsol1[i] = 
   NDSolveValue[{D[u[t, x, y], t] == 

      D[u[t, x, y], x, x] + D[u[t, x, y], y, y] + p[i][x, y], 
     u[0, x, y] == 0, u[t, x, 0] == 0, u[t, 0, y] == 0, 
     u[t, x, 1] == 0, u[t, 1, y] == 0}, 
    u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}]; 
  a[i] = (NIntegrate[
      p[i][x, y]^2, {x, 0, 1}, {y, 0, 
       1}])/(NIntegrate[(nsol1[i][1, x, y])^2, {x, 0, 1}, {y, 0, 1}]);
  g[i] = Interpolation[
    Flatten[Table[{{x, y}, g[i - 1][x, y] - a[i]*(p[i][x, y])}, {x, 0,
        1, np1}, {y, 0, 1, np1}], 1]]; npr = i; 

  If[Abs[g[i][.5, .5] - g[i - 1][.5, .5]] <= 10^-5, Break[]], {i, 1, 
   n}]]
{Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, PlotLabel -> "Exact Source"], 
 Plot3D[g[n][x, y], {x, 0, 1}, {y, 0, 1}, 
  PlotLabel -> "Constructed Source", PlotRange -> All]}

fig2


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

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more