# Problem with solving a linear ODE with a parameter

I have been struggling with the following ode for several weeks, but I fail to solve it. So I post this question and my attempts and expect someone here could give some suggestions.

First, define an operator and a function:

op[r_] = (D[#, {r, 2}] - 1/r*D[#, r]) &;
f0[r_] = 1/2*(2 r^2 Log[r/b] - r^2 + b^2);


Next, the ode and 4 boundary conditions:

eqn = op[r][op[r][f1[r]]] == -I*(4 - 2*G*(1 + b)^2*Log[r/b] + G/2*(r^2 - b^2)*(1 + (1 + b)^2/r^2));

f1[b] == 0;
f1'[b] == 0;
bc3 = f1''[r] - 1/r*f1'[r] + G*(c*f0[r] - f1[r]);
bc4 = f0'[r] - I*(f1'''[r] - 1/r*f1''[r] + 1/r^2*f1'[r]) + f0[r]/(1 + b)^2;


where both G and b are related to a constant g (0<g<1, see below), bc3 and bc4 are the left-hand sides of the actual boundary conditions both of which vanish at r=1+b (b>0, see below). Note that c in bc3 is an unknown constant (maybe complex-valued) to be determined. There are also 2 simple bcs at r=b, please see below. The ode is defined from r=b=g/(1 - g) to r=b+1=1/(1 - g).

My goal is to find c and f1[r]. At the very least, I need to find the unknown parameter c. And this is the answer for checking

c=I*G/4(G/192*(1+2*b)^2*(41+82*b+30*b^2)-1/8*(1+2*b+5*b^2)*(3+6*b+2*b^2)-3/4*(1+b)^2*(1+2*b)*Log[(1+b)/b]+(1+b)^4*(Log[(1+b)/b])^2)-I/(2*(1+b))*(1/8*((1+2*b)/(1+b))^2-1/G)


Method 1: direct method

sol1 = DSolve[{eqn /. G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1) /. b -> g/(1 - g), (f1[b] /. b -> g/(1 - g)) == 0, (f1'[b] /. b -> g/(1 - g)) == 0,
(bc3 /. r -> (1 + b) /. b -> g/(1 - g) /. G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1)) == 0,
(bc4 /. r -> (1 + b) /. b -> g/(1 - g)) == 0}, f1, (*{r, g/(1 - g), 1/(1 - g)}*) r,
Assumptions -> g > 0 && g < 1 && c != 0]


I got an error with the unsolved system returned:

DSolve::alliv: The function f1[r] was specified without dependence on all the independent variables. Each function must depend on all the independent variables.

Update on Feb 4 according to Michael's suggestion: The warning might arise from specifying a range of independent. By replacing {r, g/(1 - g), 1/(1 - g)} by r, method 1 produces a solution.

Method 2: find a general solution firstly, then use the bcs to determine the unknown constants

f1genesol = DSolve[eqn, f1, r, Assumptions -> b > 0 && {b, G} \[Element] Reals] //FullSimplify

(*{{f1 -> Function[{r}, C[4] - 1/384 I r^2 (12 b^2 G + 24 b^3 G + 12 b^4 G - 120 r^2 - 66 G r^2 - 132 b G r^2
- 51 b^2 G r^2 + G r^4 + 192 I C[1] - 96 I C[2] + 96 I r^2 C[3] + 12 ((8 + G + 2 b G) r^2 + 16 I C[2]) Log[r] +
24 b^2 (1 + b)^2 G Log[r]^2 + 60 (1 + b)^2 G r^2 Log[r/b] - 24 (1 + b)^2 G r^2 Log[r/b]^2)]}}*)


evaluating the boundary conditions with the general solution and storing in bcs

BCs = {f1[b], f1'[b], (bc3 /. r -> 1 + b), (bc4 /. r -> 1 + b)};
bcs = BCs /. f1genesol[[1]];


using the bc3 in bcs to determine the parameter c

csol = Simplify[Solve[bcs[[3]] == 0, c], b > 0]
(*{{c -> (-384 I - 768 I b - 384 I b^2 + 96 I b^2 G + 192 I b^3 G + 96 I b^4 G + 65 I G^2 + 390 I b G^2 + 948 I b^2 G^2 +
1192 I b^3 G^2 + 813 I b^4 G^2 + 282 I b^5 G^2 + 38 I b^6 G^2 + 192 G C[1] + 384 b G C[1] + 192 b^2 G C[1] - 384 C[2] -
96 G C[2] - 192 b G C[2] - 96 b^2 G C[2] - 768 C[3] - 1536 b C[3] - 768 b^2 C[3] + 96 G C[3] + 384 b G C[3] + 576 b^2 G C[3]
+ 384 b^3 G C[3] + 96 b^4 G C[3] + 384 G C[4] - 12 I (1 + b)^4 G (-16 + 5 (1 + b)^2 G) Log[1 + 1/b] +
24 I (1 + b)^4 G (-8 + (1 + b)^2 G) Log[1 + 1/b]^2 - 12 I (1 + b)^2 (-64 + (1 + b)^2 (1 + 2 b) G^2 + 16 I G C[2]) Log[1 + b] -
24 I b^2 (1 + b)^4 G^2 Log[1 + b]^2)/(192 G (-1 - 2 b + 2 (1 + b)^2 Log[1 + 1/b]))}}*)


Substituting c into bc3 and solving for the constants C[i] together with the remaining bcs

Solve[bcs[[1]] == 0 && bcs[[2]] == 0 && (bcs[[3]] /. csol[[1]]) == 0 && bcs[[4]] == 0, {C[1], C[2], C[3], C[4]}]


But, it gives null.

Method 3: find solution with first 2 bcs, then solve the remaining two unknown constants with Resolve:

sol3 = DSolve[{eqn /. G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1) /. b -> g/(1 - g), (f1[b] /. b -> g/(1 - g)) == 0, (f1'[b] /. b -> g/(1 - g)) == 0}, f1, r, Assumptions -> g > 0 && g < 1] // Simplify

(*Not show the long solution, which includes unknowns C1, C2, C3, C4*)


next, using Resolve with Exists (Note it runs more than 6 hours before manually aborting)

Resolve[Exists[{C[1], C[2]},
Simplify[(bc3 /. r -> (1 + b) /. b -> g/(1 - g) /. G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1)) /. sol3[[1]]] == 0 && Simplify[(bc4 /. r -> (1 + b) /. b -> g/(1 - g)) /. sol3[[1]]] == 0], Complexes]


Any suggestions are welcome!

• Method 1, using r instead of {r, g/(1 - g), 1/(1 - g)} that is interpreted as specifying three variable (I think), gives me a solution. Can you use it? Commented Feb 4, 2022 at 0:59
• @Enter Probably parameter c could be specified from some condition like Norm[f1]<Infinity for any 0<g<1. In this regards could you provide some link to the original problem explanation? Commented Feb 4, 2022 at 4:07
• @MichaelE2 thank you for the suggestion. I hope it could work, but it gives a null solution for c. Please see my update. Commented Feb 4, 2022 at 11:03
• @AlexTrounev I have double-checked the problem and do not see any extra stuff. This is a question that I rearranged and extracted from multiple problems which are not easy to be combined and linked here. Commented Feb 4, 2022 at 14:07
• DSolve returns a solution such that bc3 simplifies to 0; that is, it satisfies the boundary condition. Thus solving bc3 == 0 is thus equivalent to solving 0 == 0 when the solution is plugged in, and all values of c satisfy it. That is what Solve[] indicates by returning {{}}. (No solution is indicated by {}. Technically {{}} indicates the "solution set is full dimensional" according to the docs, not necessarily that all values are solutions.) Commented Feb 4, 2022 at 17:43

This problem could be clarified as follows

op[r_] = (D[#, {r, 2}] - 1/r*D[#, r]) &;
f0[r_] = 1/2*(2 r^2 Log[r/b] - r^2 + b^2);
eqn = op[r][
op[r][f1[r]]] == -I*(4 - 2*G*(1 + b)^2*Log[r/b] +
G/2*(r^2 - b^2)*(1 + (1 + b)^2/r^2));

bc1 = f1[b] == 0;
bc2 = f1'[b] == 0;
bc3 = f1''[r] - 1/r*f1'[r] + G*(c*f0[r] - f1[r]) /. {r -> b + 1};
bc4 = f0'[r] - I*(f1'''[r] - 1/r*f1''[r] + 1/r^2*f1'[r]) +
f0[r]/(1 + b)^2 /. {r -> b + 1};

sol = DSolve[{eqn, bc1, bc2, bc3 == 0, bc4 == 0}, f1[r], r];

ff = f1[r] /. sol[[1]] /. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1),
b -> g/(1 - g)} // Simplify


Now if we evaluate Denominator[ff] // FullSimplify we got 0. We also can evaluate

Denominator[ff]

Out[]= 192 (-1 + g)^4 (-1 + g^2 - 2 Log[g]) (Log[1/(1 - g)] +
Log[g] - Log[g/(1 - g)])


And we see that (Log[1/(1 - g)] + Log[g] - Log[g/(1 - g)])=0 for 0<g<1. Therefore we have a problem 1/0 on solution ff. To avoid this problem we can first define solution with 4 arbitrary constants and evaluate bc

s4 = DSolve[eqn, f1, r];

equ1 = bc1 /. s4[[1]];

equ2 = bc2 /. s4[[1]];

equ3 = bc3 /. s4[[1]];

equ4 = bc4 /. s4[[1]];


If we try to define constants using linear equations we have

{vec, mat} =
CoefficientArrays[{equ1, equ2, equ3 == 0, equ4 == 0}, {C[1], C[2],
C[3], C[4]}]


Please, note, that vec in this case has 4 nonzero elements. Therefore solution exits only if determinant not equals zero. Determinant can be computed as

det = Det[mat // Normal] // FullSimplify

Out[]= I b (1 + b) (4 + G + 2 b G +
2 (1 + b)^2 G (Log[b] - Log[1 + b]))



For required b, G we have

det /. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1),
b -> g/(1 - g)} // FullSimplify

Out[]= 0


It means that there is no unique solution with these boundary conditions. But we can use c to define unique solution as follows

ss4 =
Solve[{equ1, equ2, equ3 == 0, equ4 == 0}, {C[1], C[2], C[3], c}];

c4 = c /. ss4[[1]] /. C[4] -> 0 // FullSimplify

Out[]= -((I (48 (10 + b (48 + b (83 + 2 b (32 + b (11 + b))))) -
4 (1 + b)^2 (-9 +
b (-48 + b (-90 + b (-66 + (-11 + b) b)))) G - (1 +
b)^5 (1 + 2 b) (71 + b (142 + 61 b)) G^2 +
2 (1 + b)^2 (12 b^2 (1 + b)^3 G Log[
1 + 1/b]^2 (4 - (1 + b)^2 G +
2 (1 + b)^2 G (Log[b] - Log[1 + b])) +
6 Log[1 + 1/b] (-8 (6 + b (16 + 3 b (5 + 2 b))) -
2 (1 + b) (5 + b (21 + 2 b (17 + 2 b (6 + b)))) G + (1 +
b)^3 (5 + b (20 + b (31 + b (22 + 5 b)))) G^2 -
2 b^2 (1 + b)^2 (1 + 2 b) G (2 + G + b G) (Log[b] -
Log[1 + b])) +
b^2 G (-12 (3 + 2 b (4 + b (3 + b))) - (1 + b)^3 (36 +
b (72 + 31 b)) G +
12 (1 + b)^5 G (Log[b] - Log[1 + b])) (Log[b] -
Log[1 + b]))))/(192 (1 + b)^3 G (-1 - 2 b +
2 (1 + b)^2 Log[1 + 1/b])))


Note c4 equals to c proposed by Enter on the line {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1), b -> g/(1 - g)},C[4]->0}. Therefore, solution f1 has a form

fsol = First[f1 /. s4[[1]] /. ss4 /. C[4] -> 0 // Simplify]


We can plot this function with r and g as follows

f1g = fsol /. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1),
b -> g/(1 - g)} // Simplify;

Table[Plot[Im[f1g[[2]]], {r, g/(1 - g), g/(1 - g) + 1},
PlotLabel -> g], {g, .1, .9, .1}]


Note that evaluation (c4 - c)/. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1), b -> g/(1 - g)} // FullSimplify not gives zero, but gives

(I (Log[1/g] + Log[g]) (5 + 39 g^2 - 63 g^4 + 19 g^6 +
48 Log[1/g]^2 (-1 + g^2 - 2 Log[g]) +
12 Log[1/
g] (-5 + 3 g^2 + g^4 - 4 Log[g]) (-1 + g^2 - 2 Log[g]) +
12 Log[g] (-10 + 21 g^2 - 6 g^4 +
2 (-4 + g^2) Log[g])))/(24 (-1 + g)^2 (-1 + g^2 +
2 Log[1/g]) (-1 + g^2 - 2 Log[g])^2)


The first multiplier in this expression is (Log[1/g] + Log[g])=0 for 0<g<1, but Mathematica not evaluates it to 0 with //FullSimplify.We can add assumption to get 0 as follows

 FullSimplify[(c4 - c)/. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1), b -> g/(1 - g)}, Assumptions -> g > 0]

Out[]= 0


• is the state "...only if determinant not equals zero..." correct? The system of homogenous linear equations represented by mat has a non-trivial solution when the matrix has zero determinant Commented Feb 4, 2022 at 15:35
• @Enter Thank you for correct expression for c. System of linear equations {equ1, equ2, equ3 == 0, equ4 == 0} is not homogenous since vec in {vec, mat} = CoefficientArrays[{equ1, equ2, equ3 == 0, equ4 == 0}, {C[1], C[2], C[3], C[4]}] has 4 nonzero elements. Commented Feb 4, 2022 at 22:48
• how can we see "c4 equaling to c" explicitly, I try (c4 - c)/. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1), b -> g/(1 - g)} // FullSimplify which does not give 0 Commented Feb 5, 2022 at 5:47
• @Enter It gives long expression with (Log[1/g] + Log[g]) multiplier that equals zero for real g. Please, see update to my answer. Commented Feb 5, 2022 at 8:21
• Sir, sorry for a further question: can Mathematica simplify c4 to the simpler form of c as given? I have tried Factor and manual derivation but still can't reduce it. Thank you so much! Commented Feb 6, 2022 at 13:46

This is an extended comment rather than a full answer.

In the method 2 part of the OP you stated But, it gives null.. Ok, so let's run a quick diagnostic check.

op[r_] = (D[#, {r, 2}] - 1/r*D[#, r]) &;
f0[r_] = 1/2*(2 r^2 Log[r/b] - r^2 + b^2);

eqn = op[r][
op[r][f1[r]]] == -I*(4 - 2*G*(1 + b)^2*Log[r/b] +
G/2*(r^2 - b^2)*(1 + (1 + b)^2/r^2));
f1[b] == 0;
f1'[b] == 0;
bc3 = f1''[r] - 1/r*f1'[r] + G*(c*f0[r] - f1[r]);
bc4 = f0'[r] - I*(f1'''[r] - 1/r*f1''[r] + 1/r^2*f1'[r]) +
f0[r]/(1 + b)^2;


This is the important bit. We use the value for the constant c that you provided as the right answer in the OP.

c = I*G/4 (G/192*(1 + 2*b)^2*(41 + 82*b + 30*b^2) -
1/8*(1 + 2*b + 5*b^2)*(3 + 6*b + 2*b^2) -
3/4*(1 + b)^2*(1 + 2*b)*
Log[(1 + b)/b] + (1 + b)^4*(Log[(1 + b)/b])^2) -
I/(2*(1 + b))*(1/8*((1 + 2*b)/(1 + b))^2 - 1);


So, we want to DSolve with all boundary conditions and the answer for c to see if we can get a consistent answer.

DSolve[{eqn,
f1[b] == 0, (f1'[r] /. r -> b) == 0, (bc3 /. r -> b + 1) ==
0, (bc4 /. r -> b + 1) == 0}, f1[r], {r},
Assumptions -> b > 0 && {b, G} ∈ Reals] // Flatten //
Factor // Simplify


The above returns an answer that does not appear to have any issues. It is what I call f11[r]. Since the output is too long and it does not take to execute this in Mathematica, I do not provide its explicit expression.

Returning to the task at hand. You have a fourth order differential equation and four boundary conditions. If you attempt to do a DSolve with all you will get a “nice” answer. So, we want to execute:

DSolve[{eqn,
f1[b] == 0, (f1'[r] /. r -> b) == 0, (bc3 /. r -> b + 1) ==
0, (bc4 /. r -> b + 1) == 0}, f1[r], {r},
Assumptions -> b > 0 && {b, G} ∈ Reals] // Flatten //
Factor // Simplify;


The output of the above computation is what I call f1[r]. Please note, that I did try to use FullSimplify but it was taking too long to finish.

Let me pose my question at this point. The constant c that you want to obtain and for which you provided a value to check appears in the expression for bc3. The only known thing about bc3 is that it vanishes at r=b+1, so it does not leave much to proceed. The constant c appears, also, in the solution for f1[r], but again the only information is that the function and its derivative vanish at r=b. I am wondering if there's any extra stuff that perhaps you forgot to include in the OP and would allow us to proceed. In any case, you can have the solution for f1[r] without any undetermined coefficients.

Finally, one can try to compare the two solutions we discussed above, namely

(f1[r] // Factor // Simplify) - (f11[r] // Factor // Simplify) //
Factor // Simplify // FullSimplify


yields

(I G (-1 - 2 b +
2 (1 + b)^2 Log[
1 + 1/b]) (48 (7 + 12 b + 4 b^2 + 16 I (1 + b)^3 c) -
24 (1 + b)^3 (3 + 2 b (3 + b)) (1 + b (2 + 5 b)) G + (1 +
b)^3 (41 + 82 b + 30 b^2) (G + 2 b G)^2 +
48 (1 + b)^5 G Log[
1 + 1/b] (-3 - 6 b + 4 (1 + b)^2 Log[1 + 1/b])) ((b - r) (b +
r) + 2 r^2 (-Log[b] + Log[r])))/(1536 (1 + b)^3 (4 + G +
2 b G + 2 (1 + b)^2 G (Log[b] - Log[1 + b])))

• before executing the 2nd DSolve, we need to Clear[c] first, right? Commented Feb 4, 2022 at 13:43
• Thank you for the help! may i comment on "The only known thing about bc3 is that it vanishes at r=b+1, ... but again the only information is that the function and its derivative vanish at r=b,...". Actually, bc4 relates the 1st, 2nd and 3rd-order derivatives of f1 and thus involves c implicitly, so can we proceed with bc4 somehow? Btw, I am sorry for a typo in the last term of c (updated). Commented Feb 4, 2022 at 13:52
• I have double-checked the problem and do not find any other useful information... Well, as Alex mentioned, a physically meaningful solution should be bounded Norm[f1]<Infinity. Btw, could you suggest any kind of extra conditions that would allow us to proceed? Commented Feb 4, 2022 at 14:01
• If we evaluate last expression as %/. {G -> (4*(1 - g)^2)/(g^2 - 2 Log[g] - 1), b -> g/(1 - g)} // FullSimplify we got answer ComplexInfinity. See my answer for clarification. Commented Feb 4, 2022 at 14:12
• @DiSp0sablE_H3r0 with the typo in c being corrected, your last expression is evaluated to zero, which means c should have been determined by the given conditions. Am i right? Commented Feb 4, 2022 at 14:19