1
$\begingroup$

My code is creating an error when I try to evaluate an integral. Could I ask you to reproduce this code and see what is wrong? Thank you.

Here I simply define two variables:

Find SolVm1:

odeM1 = 1/2*sigma*Vm''[x] + ((delta*a0)/a1 - a1*x)*Vm'[x] - em*Vm[x] == -((2*x - 2*lambda*ep*(a0 + a1*x))/(em - ep));
SolVm1 = DSolveValue[odeM1, Vm[x], x] // FullSimplify

Find SolVp3:

odeP3 = 1/2*sigma*Vp''[x] + ((delta*a0)/a1 - a1*x)*Vp'[x] - ep*Vp[x] == -((-2*x + 2*em*(a0 + a1*x)* lambda)/(em - ep));

SolGeneralHomoP3 = C[2] Hypergeometric1F1[ep/(2 a1), 1/2, (a1^2 x - a0 delta)^2/(a1^3 sigma)]

odeP3 /. Vp -> (h0 + h1 # &) /. Equal -> Subtract;
eqn = CoefficientList[%, x];
solh = Solve[eqn == 0, {h0, h1}][[1]]
Simplify[odeP3 /. Vp -> (h0 + h1 # &) /. solh]
SolParticularNonHomoP3 = h0 + h1 x /. solh

SolVp3 = SolGeneralHomoP3 + SolParticularNonHomoP3 // FullSimplify

From here on is where I make the mistake:

Rename the constant in SolVp3 and transform SolVm1 and SolVp3 into functions:

Vm[x_] = SolVm1 /. Integrate -> Inactive[NIntegrate];
Vp[x_] = SolVp3 /. C[2] -> C[3];

Find the constants:

Jx[x_] = em*Vm[x] + ep*Vp[x];
Jxx[x_] = em*D[Vm[x], x] + ep*D[Vp[x], x];
Jpi[x_] = Vm[x] + Vp[x];
Jpipi[x_] = D[Vm[x], x] + D[Vp[x], x];

constants = Solve[{Jx[xhat] == 0, Jxx[xhat] == 0, Jpi[xhat] == (2*theta)/((epsilon - 1)*(1 + curlyphi))*xhat}, {C[1], C[2], C[3]}];

SolJx[x_] = Jx[x] /. constants[[1]];
SolJxx[x_] = Jxx[x] /. constants[[1]];
SolJpi[x_] = Jpi[x] /. constants[[1]];
SolJpipi[x_] = Jpipi[x] /. constants[[1]];

Specify parameters:

delta = 2/100;
epsilon = 3/2;
curlyphi = 5;
xi = 5/1000;
lambda = 1/2;
theta = lambda*(epsilon - 1)*(1 + curlyphi) + 0.05;
sigma = (20/100)^2;
muA = 5/100;
sigmaA = 10/100;
rn = delta + muA - sigmaA^2;
em = (delta*theta - Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*theta);
ep = (delta*theta + Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*theta);
a1 = (-delta*theta + Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*theta);
a0 = -((xi + rn - 1/2*sigma - sigmaA*Sqrt[sigma])/(delta + a1))*a1;
a = (-delta*theta + Sqrt[delta^2*theta^2 + 4*lambda*(epsilon - 1)^2*(1 + curlyphi)^2])/(2*lambda*(epsilon - 1)*(1 + curlyphi));

Find xhat (ultimate goal):

Solxhat = NSolve[Activate[SolJpipi[xhat]] == (2*theta)/((epsilon - 1)*(1 + curlyphi)), xhat]

The error:

enter image description here

What I don't understand: why does this error arise? I am setting x=xhat, this should not affect the variable of integration but only the limit.

I can notice something is wrong because this does not compute:

Activate[SolJx[-0.1]]
$\endgroup$

2 Answers 2

0
$\begingroup$

K[1] is a localized variable of integration.

Its occurrence is a proof, that an intermediate result contains an unresolved integral. mostly eg

f=InverseFunction[NIntegrate[expr,{K[1],0,K[2]]]

The easiest way to locate it: name and Print every line by (res=(procedure); Print[res]; res)

$\endgroup$
1
  • $\begingroup$ I wan't able to solve the problem trying to do this. When I attempt to print line by line the code crashes. I have created another post about it. $\endgroup$
    – NC520
    May 17, 2023 at 8:56
0
$\begingroup$

Maybe so:

 delta = 2/100;
 epsilon = 3/2;
 curlyphi = 5;
 xi = 5/1000;
 lambda = 1/2;
 theta = lambda*(epsilon - 1)*(1 + curlyphi) + 5/100;
 sigma = (20/100)^2;
 muA = 5/100;
 sigmaA = 10/100;
 rn = delta + muA - sigmaA^2;
 em = (delta*theta - 
 Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*
 theta);
 ep = (delta*theta + 
 Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*
 theta);
 a1 = (-delta*theta + 
 Sqrt[delta^2*theta^2 + 4*theta*(epsilon - 1)*(1 + curlyphi)])/(2*
 theta);
 a0 = -((xi + rn - 1/2*sigma - sigmaA*Sqrt[sigma])/(delta + a1))*a1;
 a = (-delta*theta + 
 Sqrt[delta^2*theta^2 + 
   4*lambda*(epsilon - 1)^2*(1 + curlyphi)^2])/(2*
 lambda*(epsilon - 1)*(1 + curlyphi));

 odeM1 = 1/2*sigma*Vm''[x] + ((delta*a0)/a1 - a1*x)*Vm'[x] - 
 em*Vm[x] == -((2*x - 2*lambda*ep*(a0 + a1*x))/(em - ep)) // 
 FullSimplify // Expand

 eq1 = B + F*x + G*x + H*Vm'[x] + A*x*Vm'[x] == A*Vm[x] + J*Vm''[x] + L*Vm''[x](*Simplify equation*)

 Sol = DSolveValue[eq1, Vm[x], x] // FullSimplify(*Simple General solution without Integrals*)

 (*1/(2 A^3 (J + L)) (A^(3/2)  Sqrt[
J + L]  (2  Sqrt[
   A]  (B Sqrt[J + L] + F Sqrt[J + L] x + G Sqrt[J + L] x + 
     Sqrt[2]  Sqrt[A]  (H + A x)  C[1] + 
     A  E^((H + A x)^2/(2 A (J + L)))  Sqrt[J + L]  C[2]) - 
  Sqrt[2 \[Pi]]  (E^((H + A x)^2/(
      2 A (J + L))) (F + G) (J + L) Erf[(H + A x)/(
       Sqrt[2] Sqrt[A] Sqrt[J + L])] + 
     A  (H + A x)  C[
      2]  Erfi[(H + A x)/(Sqrt[2] Sqrt[A] Sqrt[J + L])])) + (F + 
  G) (H + A x)^3 HypergeometricPFQ[{1, 1}, {3/2, 2}, (H + A x)^2/(
 2 A (J + L))])*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.