Solving Differential Algebraic Equations as BVP

Question

I am trying to solve a set of DAEs.

\begin{equation} -4 \nu (\lambda(s))^{(-1 - 4 \nu)} \theta'(s) \lambda'(s) + (\lambda(s))^{(-4 \nu)} \theta''(s) = -\alpha_y \cos\theta(s) + \alpha_x \sin\theta(s) \end{equation}

\begin{equation} (\lambda(s))^{(-2 \nu)} \log(\lambda(s)) = f_s (\alpha_x \cos\theta(s) + \alpha_y \sin\theta(s)) \end{equation}

\begin{equation} \theta(0) = 0 \end{equation}

\begin{equation} \theta'(1) = \beta \end{equation}

where $\lambda$ and $\theta$ are two variables, varying over the range $s \in [0,1]$. $\alpha_x, \alpha_y, \beta, f_s, \nu$ are constants. When I try to solve them numerically using NDSolve, Mathematica gives me an error saying DAEs must be given as IVPs.

The code I use is given below

i[s] = (lambda[s])^(-4 nu)
i'[s] = D[i[s],s]
Eqn1 = theta''[s] i[s] + theta'[s] i'[s] == alphax Sin[theta[s]] - alphay Cos[theta[s]]
Eqn2 = (lambda[s])^(-2 nu) Log[lambda[s]] == fs*(alphax Cos[theta[s]] + alphay Sin[theta[s]])
BC1 = theta[0] == 0
BC2 = theta'[1] == beta

param = {alphax->0.1, alphay->0.1, beta->0.1, nu->0.3, fs->10^-6}
thetaSol = NDSolve[{Eqn1,Eqn2,BC1,BC2}/.param,{theta,lambda},{s,0,1}]

If I can solve the second equation to obtain $\lambda(s)$ as a function of $\theta(s)$, then I can eliminate the second equation and solve it as a second order ODE in $\theta(s)$. However, I believe this sort of equation has a solution using the Lambert W function.

Can I use Mathematica to solve this system of equations?

Answer 1

NDSolve wants a IVP, why not give it one?:

bc3 = theta'[0] == slope;
thetaSol = 
 ParametricNDSolve[{Eqn1, Eqn2, BC1, bc3} /. param, {theta, lambda}, {s, 0, 1}, slope];

Then what we need to do is just looking for a proper slope that satisfies BC2, which is usually a task for FindRoot, let's have a look at the slope - theta'[1] relation as insurance:

Plot[First@BC2 /. theta -> theta[slope] /. thetaSol /. param // Evaluate, {slope, -1, 1}]

OK, theta'[1] == 0.1 is near 0, so:

slopeneeded= FindRoot[BC2 /. theta->theta[slope] /. thetaSol /. param//Evaluate, {slope, 0}]

{slope -> 0.191621}

{θ[s_], λ[s_]} = {theta[slope][s], lambda[slope][s]} /. thetaSol /. slopeneeded;

Plot[{θ[s], λ[s]}, {s, 0, 1}]

Answer 2

With your set of parameters your second equation looks as follows:

param = {alphax -> 0.1, alphay -> 0.1, beta -> 0.1, nu -> 0.3, fs -> 10^-6};
eqn2 = (lambda[s])^(-2 nu) Log[lambda[s]] == 
    fs*(alphax Cos[theta[s]] + alphay Sin[theta[s]]) /. 
   param /. {lambda[s] -> λ, theta[s] -> θ}

(*  Log[λ]/λ^0.6 == (0.1 Cos[θ] + 0.1 Sin[θ])/1000000   *)

One can see that the right-hand side has the factor of 10^-7. In principle I would say that with a perfect accuracy you can assume the right-hand side to be equal to zero. Then lambda=1, and you can easily solve your equation Nr. 1:

 thetaSol = 
 NDSolve[{theta''[s] == -0.1` Cos[theta[s]] + 
      0.1` Sin[theta[s]], theta[0] == 0, theta'[0] == 0.1}, 
   theta[s], {s, 0, 1}][[1, 1]]

Plot[Evaluate[theta[s] /. thetaSol], {s, 0, 1}]

This should return this plot:

If, however, for some reason (that I am unable to imagine) you need these tiny differences, then you may process as follows. This

lst = Table[{θ, 
 FindRoot[
   Log[λ]/λ^0.6` == (
    0.1` Cos[θ] + 0.1` Sin[θ])/1000000, {λ, 
    0.9}, WorkingPrecision -> 8][[1, 2]]}, {θ, 0, 2 π, 
 0.1}] // Quiet;

gives you a table of solutions with the structure {theta, lambda}. This fits it:

lst1 = lst /. {θ_, λ_} -> {θ, λ - 1};
ft = Fit[lst1, {1, Sin[θ], Cos[θ]}, θ]
Show[{
  lstPl = ListPlot[lst1],
  Plot[ft, {θ, 0, 2 π}]
  }]

This is the returned visualization of the fitting quality:

Note that I fit {theta, lambda-1} instead of the initial list. Otherwise the differences would be difficult to recognize. Now your result for lambda[theta]is 1+ft. That is

λ =  1.000000000 + -1.8133184287449915`*^-8 + 
     1.0016280409439516`*^-7 Cos[θ] + 9.982597207356538`*^-8 Sin[θ]

Have fun!