2
$\begingroup$

The differential equation in question.

I have a problem when I want to solve this integro-differential equation where the $\chi' (U)$ appears in the integrand. If statement gives the $chiIntegral=0$ for $umax=20$ since the upper bound of the integral becomes smaller than the lower bound so it properly solves a homogeneous differential equation.

However when i set $umax$ to higher values $10^8-10^{10}$, I get error, because then if statement tries to calculate the integral which does not work due to the following message:

NIntegrate::inumr: The integrand (chiPrime[U] SphericalBesselJ[2,1126.62 -U])/(1126.62 -U)^2 has evaluated to non-numerical values for all sampling points in the region with boundaries {{1125.74,1126.62}}.

My humble guess is that it does not understand maybe how to handle $\chi'$ inside the integral since knowing $\chi'$ requires the priori of $\chi$

Here is the code:

(Initial conditions and plot range)

umin = 10^-15; umax = 20; Clear[chiPrime];

(Initial condition is chi(umin)=1 and chi'(umin)=0)

chiInitCond = {chi[umin] == 1, chi'[umin] == 0};

(Define epsilon to have no singularity when U -> u in chiIntegral)

epsilon = 10^-25;

(Kernel:j2/(u-U)^2)

j2[x_] := SphericalBesselJ[2, x]; kernel[u_, U_] := j2[u - U]/(u - U)^2;

(Memoized integral function using chiPrime)

Clear[chiIntegral] chiIntegral[u_?NumericQ] := If[u > uNuDec + epsilon, NIntegrate[kernel[u, U]*chiPrime[U], {U, uNuDec, u - epsilon}, Method -> "GaussKronrodRule", MaxRecursion -> 20, AccuracyGoal -> 6, PrecisionGoal -> 6], 0]

(Define the first-order system)

chiEqSystem = {chi'[u] == chiPrime[u], chiPrime'[u] == -2 aURatioscaled[u] chiPrime[u] - chi[u] - 24 fNuScaled[u]*(aURatioscaled[u]^2)*chiIntegral[u]};

(Initial conditions)

chiInitCond = {chi[umin] == 1, chiPrime[umin] == 0};

(Solve)

chiSolutionNonHomog = NDSolve[Join[chiEqSystem, chiInitCond], {chi, chiPrime}, {u, umin, umax}, Method -> "StiffnessSwitching", MaxSteps -> Infinity, AccuracyGoal -> 8, PrecisionGoal -> 8];

My question is: how to fix this issue?

Improve this question
$\endgroup$
7
  • $\begingroup$ What are: aURatioscaled, fNuScaled, uNuDec ? $\endgroup$ Commented Apr 25 at 8:40
  • $\begingroup$ uNuDec is just a number, fNuScaled = 0.4/(1+aURatioscaled*10^4) while aURatioscaled is a solution of another differential equation $\endgroup$ Commented Apr 25 at 9:20
  • $\begingroup$ @noone235711 In your LaTeX formula there is a quotient a'[u]/a ! Are you using the same symbol for function and parameter? $\endgroup$ Commented Apr 25 at 10:04
  • 2
    $\begingroup$ @noone235711 Please provide runnable code! $\endgroup$ Commented Apr 25 at 10:53
  • $\begingroup$ @noone235711 Is there still interest in answering the question? $\endgroup$ Commented Apr 25 at 15:08

1 Answer 1

Reset to default
7
$\begingroup$

It's a bit of a shame that so many comments have been left without feedback. But the integro-differential equation is to "beautiful" to leave question unanswered.

Problem can be solved using galerkin-method and gaussian-quadratur(partly symbolic!):

Several system parameters still aren't known,that's why I assume some values

parameters

umin = 0; (* 10^-15 ganz schön klein... *) umax = 50 ; nn = 150; a = Function[u, 1/3]; (* quotient a'[u]/a *) fnue = Function[u, -24 10]; j2[x_] := SphericalBesselJ[2, x]; ui = Subdivide[umin, umax, nn]; (*discretization *) chii = Table[\[Chi][i], {i, 1, Length[ui]}]; (* unknowns *)

FEM (element functions)

Needs["NDSolve`FEM`"] netz = ToElementMesh[Map[{#} &, ui]] \[Phi] =Map[ElementMeshInterpolation[netz, #] &,IdentityMatrix[Length[ui]]]; \[Phi]x = Map[Derivative[1 ][#] &, \[Phi]i] ;

Gaussian quadrature(points&weights)

Get["NumericalDifferentialEquationAnalysis`"]; xwGauss = (* in jeder Teildiskretisierung wird 3-Punkt Gauss(P5[x]) verwendet *) Map[GaussianQuadratureWeights[3(*5*), #[[1]], #[[2]]] & ,Partition[ui, 2, 1]] ;

integral (right hand side)

(* Gaussquadratur zwint[u] mit 0<U<u ,Element[u,ui]*) zwint = Join[{0},Flatten@Table[ (* ui[[i]] obere Integrationsgrenze *) ipU =Function[U,Evaluate[Through[\[Phi]x[U]] . chii j2[ui[[i]] - U]/(ui[[i]] - U)^2]]; Total@Map[ipU[#[[1]]] #[[2]] &, Flatten[xwGauss[[1 ;; i]], 1]] , {i, 1, Length[ui] - 1}] // Chop]; (* integral[u]=Through[\[Phi][u]].zwint *) `zwint` is symbolic, depends on unknowns `chii`!

galerkin method

fu = Function[u,fnue[u] a[u]^2 Through[\[Phi][u]] Through[\[Phi][u]] . zwint ] intU = Total[Map[fu[#[[1]]] #[[2]] &, Flatten[xwGauss, 1]]] //Simplify // Chop; (* Term \[Chi][u]*) fu = Function[u, Through[\[Phi][u]] (Through[\[Phi][u]] . chii) ]; int\[Chi] =Total[Map[fu[#[[1]]] #[[2]] &, Flatten[xwGauss, 1]]] // Simplify// Chop; (* Term \[Chi]''[u] REM Green =green-\[Phi]x \[Phi]x.chii*) green = Through[\[Phi][umax]] Through[\[Phi]x[umax]] . chii // Chop; fu = Function[u, Through[\[Phi]x[u]] (Through[\[Phi]x[u]] . chii) ]; int\[Chi]xx =green - Total[Map[fu[#[[1]]] #[[2]] &, Flatten[xwGauss, 1]]] //Simplify // Chop; (* Term \[Chi]'[u]*) fu = Function[u, a[u] Through[\[Phi][u]] (Through[\[Phi]x[u]] . chii) ]; int\[Chi]x =2 Total[Map[fu[#[[1]]] #[[2]] &, Flatten[xwGauss, 1]]] //Simplify // Chop;

equation system [chii] & solution

eqn = int\[Chi]xx + int\[Chi]x + int\[Chi] - intU; ics = {Through[\[Phi][umin]] . chii == 1,Through[\[Phi]x[umin]] . chii == 0} // Chop; (* initial conditions*) mini = NMinimize[{eqn . eqn, ics}, chii]; Plot[ElementMeshInterpolation[netz, Values[mini[[2]]]][u], {u, umin,umax}, PlotRange -> All, AxesLabel -> {u, "\[Chi][u]"}]

enter image description here

Hope it helps!

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.