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?

  • $\begingroup$ Hi ! Please provide any relevant code, so someone can start diggin' right away, not wasting time trying to type in the equations :) $\endgroup$
    – Sektor
    Commented Oct 22, 2014 at 15:53
  • $\begingroup$ I have added the code I am currently using. Thanks! $\endgroup$
    – Venkat K V
    Commented Oct 22, 2014 at 16:44
  • $\begingroup$ The numerical values of the parameters $alpha_x, alpha_y$ etc are just as an example for using NDSolve. I need to able to solve the equations for other values too. $\endgroup$
    – Venkat K V
    Commented Oct 23, 2014 at 18:04

2 Answers 2


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}]

enter image description here

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}]

enter image description here

  • $\begingroup$ I am not able to get the result you get when you solve for slopeneeded. Can you please verify it? $\endgroup$
    – Venkat K V
    Commented Oct 23, 2014 at 19:23
  • $\begingroup$ @VenkatKV Of course: i.sstatic.net/wPNna.jpg What do you mean by "ot able to get the result"? The result is wrong, or it doesn't give any result, or something else? BTW, ParametricNDSolve is a v9 function. $\endgroup$
    – xzczd
    Commented Oct 24, 2014 at 2:57
  • $\begingroup$ This is the code I am using. Similar to yours. You can see the results I get. imgur.com/ohwXcUN I notice the you are using theta'[0]=slope instead of == in your code. Maybe that's an error? $\endgroup$
    – Venkat K V
    Commented Oct 24, 2014 at 15:57
  • $\begingroup$ @VenkatKV No, it's just the StandardForm of ==. (Notice its length. ) Which version of Mathematica are you using? What's your platform? I'm using v9.0.1, vista 32bit. You can try narrowing down the scope of FindRoot, for example {slope, 0.1, 0.3} etc. $\endgroup$
    – xzczd
    Commented Oct 25, 2014 at 4:06
  • $\begingroup$ I am using version 10.0.1. Providing the range for FindRoot seems to have done the trick for the calculation of slope. I will try doing the calculations for a large range of loads to see if the method works. Thanks! $\endgroup$
    – Venkat K V
    Commented Oct 27, 2014 at 19:14

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: enter image description here

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[{θ, 
   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[θ]}, θ]
  lstPl = ListPlot[lst1],
  Plot[ft, {θ, 0, 2 π}]

This is the returned visualization of the fitting quality: enter image description here

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!

  • $\begingroup$ I think you may be interested in this site :) $\endgroup$
    – xzczd
    Commented Oct 23, 2014 at 10:10
  • $\begingroup$ Thanks for the reply, but the numerical values of $alphax, alphay$ etc are just as an example. I need to be able to solve the equations for a much larger range of values, in which case $lambda$ is not equal to $1$. However, the second part of your reply on how to numerically obtain $lambda(theta)$ is very useful. $\endgroup$
    – Venkat K V
    Commented Oct 23, 2014 at 18:00

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.