# Solving Differential Algebraic Equations as BVP

I am trying to solve a set of DAEs.

$$-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)$$

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

$$\theta(0) = 0$$

$$\theta'(1) = \beta$$

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?

• Hi ! Please provide any relevant code, so someone can start diggin' right away, not wasting time trying to type in the equations :) Commented Oct 22, 2014 at 15:53
• I have added the code I am currently using. Thanks! Commented Oct 22, 2014 at 16:44
• 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. Commented Oct 23, 2014 at 18:04

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


• I am not able to get the result you get when you solve for slopeneeded. Can you please verify it? Commented Oct 23, 2014 at 19:23
• @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. Commented Oct 24, 2014 at 2:57
• 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? Commented Oct 24, 2014 at 15:57
• @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. Commented Oct 25, 2014 at 4:06
• 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! 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:

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!

• I think you may be interested in this site :) Commented Oct 23, 2014 at 10:10
• 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. Commented Oct 23, 2014 at 18:00