1
$\begingroup$

Is it normal that Mathematica takes more than 2 hours to solve this system of differential equations?

ODEs:

OdeM = 1/2*\[Sigma]X^2*V1''[x] + ((a0*\[Delta])/a1 - a1*x)*V1'[x] - em*V1[x] == -((2*x - 2*ep*(a0 + a1*x)* \[Lambda])/(em - ep));
OdeP = 1/2*\[Sigma]X^2*V2''[x] + ((a0*\[Delta])/a1 - a1*x)*V2'[x] - ep*V2[x] == -((-2*x + 2*em*(a0 + a1*x)* \[Lambda])/(em - ep));

Initial conditions:

Ic1 = em*V1[xhat] + ep*V2[xhat] == 0;
Ic2 = em*V1'[xhat] + ep*V2'[xhat] == 0;
Ic3 = V1[xhat] + V2[xhat] == (2*\[Theta])/((\[Epsilon] - 1) (1 + \[CurlyPhi]))*xhat;
Ic4 = V1'[xhat] + V2'[xhat] == (2*\[Theta])/((\[Epsilon] - 1) (1 + \[CurlyPhi]));

Solver:

Sol = DSolveValue[{OdeM, OdeP, Ic1, Ic2, Ic3, Ic4}, {V1[x], V2[x]}, x] //FullSimplify

The same system without the four initial conditions takes barely any time. I am afraid I am doing something wrong, even though it seems a pretty standard problem.

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
3
$\begingroup$

The same system without the four initial conditions takes barely any time.

It has nothing to do with IC. It is FullSimplify which takes very long time. Try

Timing[DSolveValue[{OdeM, OdeP, Ic1, Ic2, Ic3, Ic4}, {V1[x], V2[x]},x]]

On my PC, V 13.2.1 it takes 4.5 seconds

Mathematica graphics

ps. I renamed your σX^2 to $w$. As I mentioned in your other question wrong-sign-in-variation-of-parameters-method , it is bad idea to use two letters for one variable. Why do you insist in using σX for a variable instead of say w or any other single letter?

Any way, it is well known that FullSimplify can slow things alot as it tries many types of simplifications. Try to use it only if needed. Simplify is much faster than FullSimplify

Full code

OdeM = 1/2*w^2*V1''[x] + ((a0*δ)/a1 - a1*x)*V1'[x] - 
    em*V1[x] == -((2*x - 2*ep*(a0 + a1*x)*λ)/(em - ep));
OdeP = 1/2*w^2*V2''[x] + ((a0*δ)/a1 - a1*x)*V2'[x] - 
    ep*V2[x] == -((-2*x + 2*em*(a0 + a1*x)*λ)/(em - ep));

Ic1 = em*V1[xhat] + ep*V2[xhat] == 0;
Ic2 = em*V1'[xhat] + ep*V2'[xhat] == 0;
Ic3 = V1[xhat] + V2[xhat] == (2*θ)/((ϵ - 1) (1 + φ))*xhat;
Ic4 = V1'[xhat] + V2'[xhat] == (2*θ)/((ϵ - 1) (1 + φ));
Timing[DSolveValue[{OdeM, OdeP, Ic1, Ic2, Ic3, Ic4}, {V1[x], V2[x]}, x]]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you. I used FullSimplify because I was hoping to find in a clear way what the four constants are. $\endgroup$
    – NC520
    Commented Apr 22, 2023 at 16:56
3
$\begingroup$ 
gensol = DSolve[{odeM, odeP(*,ic1,ic2,ic3,ic4*)}, {V1, V2}, x];

linsys = 
 CoefficientArrays[{ic1, ic2, ic3, ic4} /. Equal -> Subtract /. 
   First@gensol, Array[C, 4]]
paramsol = Thread[Array[C, 4] -> LinearSolve @@ Reverse[{-1, 1} linsys]];
icsol = gensol /. paramsol

Mathematica graphics

{ic1, ic2, ic3, ic4} /. icsol // Simplify

{{True, True, True, True}}  *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ This is not running for me $\endgroup$
    – NC520
    Commented Apr 22, 2023 at 16:55
  • $\begingroup$ @NC520 I changed the initial upper case to lower case (following best practices). Except I forgot to do it for V1 and V2. Or: I forgot to change a sign (now fiexed). Otherwise, you'll have to be clearer. It works for me. $\endgroup$
    – Michael E2
    Commented Apr 22, 2023 at 17:05
  • $\begingroup$ It works now. Thank you. $\endgroup$
    – NC520
    Commented Apr 22, 2023 at 17:08

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.