5
$\begingroup$

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?

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

2 Answers 2

Reset to default
4
$\begingroup$

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

Improve this answer
$\endgroup$
7
  • $\begingroup$ I am not able to get the result you get when you solve for slopeneeded. Can you please verify it? $\endgroup$ 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$ 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$ 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$ 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$ Commented Oct 27, 2014 at 19:14
1
$\begingroup$

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[{θ, 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: 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!

Improve this answer
$\endgroup$
2
  • $\begingroup$ I think you may be interested in this site :) $\endgroup$ 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$ Commented Oct 23, 2014 at 18:00

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.