I have a question about a first order non-linear differential equation. I have tried many method to solve this problem but not successful yet. Here is my question;

$$w^4 + (w')^2 = g(t)$$

$$w' = \left[g(t)- w^4\right]^{0.5}$$

I have tried runge kutta 4th order and 4-5 adaptive step size and cash-carp methods up to now, I think they are very good methods to solve this IVP but still not give expected values for some places.

Do you have another suggestions to solve this? Thanks in advance

  • $\begingroup$ Do you have any information on $g$? Does it fit in this direction? From the associated problem, do you know that there must be a real valued solution? $\endgroup$ – Nikolaj-K Feb 20 '13 at 14:43
  • 2
    $\begingroup$ Pleas use TeX-like formatting for math. For more information on this type of formatting, read the FAQ. $\endgroup$ – user41994 Feb 20 '13 at 14:46
  • 2
    $\begingroup$ From the two comments, you seem to have quite a bit more information that you provided in the initial question. I feel typing it out can't be a bad idea. $\endgroup$ – Nikolaj-K Feb 20 '13 at 14:57
  • 1
    $\begingroup$ Jacob, if $g$ falls abruptly, then the solution may cease to exist. Moreover, the numerical methods become unstable even before that. You may want to look at artofproblemsolving.com/Forum/viewtopic.php?f=73&t=216577 where my way out was not to solve it at all. Now, if you have full $g$ data (so you can look well ahead), there are better ways, but if you do it in real time, the difficulty is real, not just technical. $\endgroup$ – fedja Feb 20 '13 at 15:09
  • 1
    $\begingroup$ Yes, when you have a decent function that is far away from the breakdown regimes, the standard methods work. The problem there was that once the centrifuge was running, you had to do something even when the theory said that there was no solution, so you should either restrict the input and use high accuracy methods on the descents, or just "do something different" when the behavior gets close to critical. Actually, I'm happy to return to this question, so if you have some problem with the centrifuge equation that cannot be solved by mere increase in precision, let me know :). $\endgroup$ – fedja Feb 20 '13 at 17:33

OK, here are some (unpleasant, but apparently correct) observations. We are still trying to solve $w^4+\dot w^2=H$ in real time. I'll also assume that you want $w\ge 0$, i.e., you do not want to change the direction of your rotation in the process. Note that allowing the direction change won't help much in the sense that the final conclusion will stay the same but the argument will get more involved.

Recall (see my AoPS post) that invariant quantities are $\dot H/H^{5/4}$ and $\ddot H/H^{3/2}$. Assume that we want to create a control satisfying $w^4+\dot w^2\in [(1-\delta)H,(1+\delta)H]$ for all reasonable regimes $H$. Suppose that the pilot is not allowed to change $H$ or $\dot H$ abruptly but he can change $\ddot H$ any time and in any way within some reasonable limits, so $\ddot H$ does not need to be continuous. I believe, this is more or less what happens in reality but if I am wrong here, correct me.

Observation 1. Ascending regimes are asymptotically stable. What it means is that if $H$ is any ascending regime that starts at some distant past time moment $-T$ and equals some constant $H_0$ up to moment $0$ (we shall call such regimes canonical), then any non-negative solution of the equation $w^4+\dot w^2=H$ that is defined on the entire interval $[-T,T]$ has essentially the same values $w(t),\dot w(t)$ as the canonical solution whose initial condition is given by $w=H_0^{1/4}$ on $[-T,0]$.

Observation 2. Suppose that $H$ is a canonical ascending regime and $w$ is the corresponding canonical solution with $\dot w\ge 0$. Suppose also that $w_1$ satisfies $w_1^4+\dot w_1^2<H$ on $[-T,T]$. Then $w_1\le w$. Indeed, we clearly have $w_1\le w$ on $[-T,0]$. On the other hand, the graphs can never cross because at the first moment when $w_1=w$, we must have $\dot w_1<\dot w$, so a tiny bit before that we must have $w_1>w$, which contradicts the assumption that our crossing moment is the first one.

Observation 3. Suppose that $H$ is a canonical ascending regime and $w$ is the corresponding canonical solution with $\dot w\ge 0$. Suppose also that $c>0$ (it can be arbitrarily small) and $w_1$ satisfies $w_1^4+\dot w_1^2>(1+c)H$ with on $[-T,T]$. Then $w_1\ge w$ on $[0,T]$ if $T>T(c)$. Indeed, we have $w_1(0)\approx (1+c)^{1/4}H_0^{1/4}>w(0)$. Again, consider the first crossing moment if it exists. At that moment we must have $|\dot w_1|>\dot w$. If $\dot w_1>0$, we get the contradiction in the same way as before. However, if $\dot w_1<0$, $w_1$ must continue to head down and its derivative can only grow in absolute value in that process until the direction gets reversed so this situation corresponds to a prohibited regime.

Observation 4. The time in the equation is reversible, so instead of investigating the descending regimes, we can investigate the ascending ones. The question we'll ask will be the following. Take any canonical regime $H$. Suppose we have a clever control $v$ such that $v(t)$ depends only on $H(s)$ with $s\ge t$ (there are no assumptions about how you designed it: solving the ODE, writing some explicit formula, asking an oracle, whatever). All we know is that our control doesn't create relative error more than $\delta$ for every regime anywhere and it is not allowed to see the past (remember that the time is now running backwards). What can we say about the value $v(t)$ for some fixed $t>0$?

The catch is that now we can glue any canonical regime on $[-T,t]$ which has correct values $H(t)$ and $\dot H(t)$ to the piece we have on $[t,T]$ and $v(t)$ is determined by that piece alone. Thus, $v(t)$ should make sense for any piece we glue.

Observation 5. Let us glue a canonical piece $H=w^4+\dot w^2$ with $\dot w\ge 0$ on $[-T,t]$ (note that we can choose $w$ instead of choosing $H$). Then the solution $(1+\delta)^{1/2}w$ will give us a canonical regime that is at least $(1+\delta)H$. Similarly, $(1-\delta)^{1/2}w$ will give us a canonical regime that is at most $(1-\delta)H$. Thus, our clever control must satisfy $(1-\delta)^{1/2}w(t)\le v(t)\le (1+\delta)^{1/2} w(t)$. Note that it should hold for every admissible $w$. Now you, probably, smell trouble: we have many different $w$ to try and, if they don't give the same value, the exact solution is impossible. Moreover, if those values may be noticeably different, $\delta$ cannot be too small.

Observation 6. We still have two conditions to satisfy at $t$, which, in terms of $w$, can be written as $w^4+\dot w^2=H$, $2w^3\dot w(2+\psi)=\dot H$ where $\psi=\ddot w/w^3$. Note that we have 3 free parameters $w,\dot w,\ddot w$ and only 2 equations. Now it remains to find the quantitative estimates for the freedom we have. To this end, we will compare the scale invariant quantities $I=w^4/H$ for two regimes with the same $\dot H/H^{5/4}$. Since the controller I gave creates dismal results near the bottom of the cosine wave $2+\cos kt$ for $k>2$, it will be a good idea to look at the "quadratic departures" from the ground state and compute the maximal uncertainty along the way. For technical reasons, it'll be convenient to look at $w(t)=1+at^2$, i.e., $H(t)=(1+at^2)^4+4a^2t^2$. The value $k=2$ corresponds to $a\approx 0.35$ here. This regime gives the invariant derivative quantity $\frac{8at(1+at^2)^3+a}{[(1+at^2)^4+4a^2t^2]^{5/4}}$. This is to be compared with the time/scale invariant regime $w(t)=b/t$, i.e., $H(t)=\frac{b^4+b^2}{t^4}$ with the invariant derivative quantity $\frac{4}{(b^4+b^2)^{1/4}}$. The ratio of the corresponding $I$-values is $$ \frac{(1+at^2)^4(b^4+b^2)}{[(1+at^2)^4+4a^2t^2]b^4} $$ where $b$ is determined from the equation $$ b^4+b^2=\frac{[(1+at^2)^4+4a^2t^2]^{5}}{[8at(1+at^2)^3+a]^4} $$ I will leave it to you to run a simple script maximizing over $t$ ($t=0.2$ already gives you something to be unhappy with for $a\in(0,1)$) and to recast the corresponding "inevitable uncertainty ratio" into the $\delta$-values. The table is like that:

$a=0.1, \delta=0.006$

$a=0.2, \delta=0.018$

$a=0.3, \delta=0.036$

$a=0.4, \delta=0.058$

$a=0.5, \delta=0.085$

$a=0.6, \delta=0.116$

$a=0.7, \delta=0.150$

$a=0.8, \delta=0.188$

$a=0.9, \delta=0.228$

As you can see, at the "breakdown point" $a=0.35$, the error of $4\%$ in the $H$-profile approximation is just inevitable and my controller is about that accurate. If $a$ gets to $0.6$, which corresponds to $k=\sqrt{8(a+a^2)}=2.8$ in the cosine model, you are going to be $10\%$ off somewhere no matter what you do. So, the version of the problem you are trying to solve is not just hard but unsolvable. You have to restrict the input to something reasonable and be happy with approximate control. However, you can try to investigate what features of the current control you dislike and design a better control in some respects. Just don't waste your time trying to achieve the impossible :).

  • $\begingroup$ P.S. I suggest you show this all to ETCKerry as well. As usual, you both are welcome to ask as many questions as you want :). $\endgroup$ – fedja Feb 24 '13 at 17:29
  • $\begingroup$ Thank you fedja for your valuable time and effort again. I will compare your observations with the metdod I have tried to solve the equation. I would send the g regime data but after your observations I have decided to combine my data with your observations. With regard to my next step, I will try how I can reduce the error for omega and g as much as possible. Thanks again. You are right imposible is imposible but once upon a time it was impossible to find electricity :) $\endgroup$ – yakupc Feb 25 '13 at 9:28
  • $\begingroup$ As it became clear now, the first step in reducing the error as much as possible is to understand clearly what exactly are the restrictions on the $H$ (or, if it is easier for you to understand, on the $G$) profiles. The point is that there cannot be any good or even decent controller for harsh regimes in principle but, once you restrict harshness and tell the objectives in a precise way, we can try to find the optimal controller for those limits. Send the data anyway: I'd like to look at what you have there no matter what :). $\endgroup$ – fedja Feb 25 '13 at 12:12
  • $\begingroup$ Ok, I'll provide what I have soon.Thanks $\endgroup$ – yakupc Feb 25 '13 at 12:28


I will try to describe my way to solve this ODE accurately as much as I can.

Below is the g equation we need. To recall the problem, I will describe the variants again. r represents the radius of circular motion, g represents the gravity and they are constant.

$$G^2 = ((w^2)^*r/g)^2 + (dw/dt * r/g)^2 + 1 $$

After substituting the w from Runge Kutta 4th order for the w in the above equation to get desired g, as you know, I have not been successful yet for decreasing commands.

After many attempts, I thought the problem is tangential part of the g equation, and therefore I removed dt from the above equation(new equation is shown below), and then, I compared the exact g values with the g values from below equation, result is satisfactory, and even beyond that.

$$G^2 = ((w^2)^*r/g)^2 + (dw * r/g)^2 + 1 $$

After I get satisfactory result from above equation, I decide to take account the error in w.

$x = w(t) + error$

$dx = w(t) +error - w(t-1)$

$$(x^2*r/g)^2 + (dx/dt * r/g)^2 + 1 = (w^2*r/g)^2 + (dw * r/g)^2 + 1 $$ $$ x^4 + (dx/dt)^2 = w^4 + (dw )^2 $$

After solving this equation, reach the error equation in quartic. As you can imagine after some successful result, complex roots interrupt this dream.

Is there any way to interpret the complex roots? I asked this because, I get many complex roots as shown below. $$ r1 = a+ib$$ $$r2 = a-ib $$ $$ error(t) = exp(a*t)[C1*cos(b*t)+C2*sin(b*t)] $$

C1 and C2 are the arbitrary constants and unfortunately I couldn't select a good value for C1 and C2.

Do you think this way is deserved to continue or I will turn another way to solve?


  • 1
    $\begingroup$ As you could see from my posts on AoPS, ETCKerry had exactly the same problem with the complex roots :). I'll look at what you are doing and will be happy to think of the problem again myself (but I have to warn you that it may be a matter of weeks, not hours, and there is no warranty on the outcome). To start with, tell me a few explicit G-regimes where your approach breaks. I want to look at two things: 1) whether the real-time solution exists for them in principle and 2) whether I can get it without running into trouble. $\endgroup$ – fedja Feb 21 '13 at 14:53
  • $\begingroup$ Sorry I didn't mention it; during reviewing your post on AoPS, I realized that ETCKerry who is, and I realized that DE is same in the post you sent:) we are working same company but different location :) $\endgroup$ – yakupc Feb 21 '13 at 16:35
  • $\begingroup$ I've checked DE you prepared,It is good for low decreasing g command. However It is not working for high G decreasing commands. while I was observing the G'(derivative of the resultant g), the data is broken end of the function;with regard to this problem I will send some data later.Anyway I am so interest in math/DE and love it,therefore I have started to try to solve this problem for high g decreasing commands again myself.I'll also provide data where g breaks by using RK4 and RK CashCarp adaptive step size methods. Thank you in advance for your help.I'll call Kerry this conversation again. $\endgroup$ – yakupc Feb 21 '13 at 16:38
  • 1
    $\begingroup$ Great! I'll wait for the data on which both your and my methods break. My own favorite "supercritical" test regime is $H=2+\cos(kt)$. When $k$ gets above $2$, my controller outputs a terrible result near the bottom of the $\cos$ wave but I suspect that you may have the same trouble with your methods here. As I said, I'm happy to return to the problem. Just give me some test regimes that are still problematic and remember that some regimes (a very quick offset from the constant, say) are just impossible no matter what method you use and all we can hope for with them is "something reasonable". $\endgroup$ – fedja Feb 22 '13 at 5:48

Below is the result for ODE that I have used with Runge Kutta Cash-Karp(RKCK) with adaptive step size method.

Blue is the $desiredG$, green is the $resultG$, red is the $desiredG'$ and black one is the $resultG'$.

As you can see, It is very good for onsets, but we can't say same thing for the offset.

method 1 http://ycengiz.files.wordpress.com/2013/02/method-1.png

And, here is the result for the ODE you have prepared. color meanings are same as above. (note that $dt$ is $0.1$, It effects the result as far as I experienced)

method 2 http://ycengiz.files.wordpress.com/2013/02/method-2.png

Eventually, both results seem that almost same for the g offsets. However onset result is better when I have taken the error(t) into account in RKCK method.

Can you explain again why we need to use $\sqrt{dt}$ instead of $dt$ or different step. When I set $dt$ to $0.01$, $g$ changes in time turn to noisy data like shown below.

method 2 with noise http://ycengiz.files.wordpress.com/2013/02/method-2-noise.png

The last chart is regarding when I changed the exact G solution. When tangential part is multiplied with $dt^2$ and compare the results I got below result. Very good offset data! It seems that DE can solve another equation (wrong equation but at least we can use it as reference function in the DE) as we expected. How can we compare these two equation and eliminate the error on omega? I will explain how error function is created myself.

method 2 with wrong solution http://ycengiz.files.wordpress.com/2013/02/method-1-with-wrong-solution.png

error(t) equation is prepared as shown below.

  • $X =$ desired $W$
  • $W = $ calculated from RK method
  • $W_{\rm old} =$ omega at previous segment
  • $dt =$ step size

$$ (X^2 \textrm{radius}/\textrm{gravity})^2 + ((X-W_{\rm old}/dt)*\textrm{radius}/\textrm{gravity})^2 + 1 \\ = (W^2\textrm{radius}/\textrm{gravity})^2 + ((W-W_{\rm old}/dt)*\textrm{radius}/\textrm{gravity})^2 dt^2 + 1 $$

we can eliminate radius and gravity from the equation, new equation is below

$$ X^4 + ((X-W_{\rm old}/dt)^2 = (W^4 + (W-W_{\rm old})^2 + 1 $$

It gives again and again complex roots same as solving G equation :)

Here is another result when I multiplied omega with a constant from error(t), I think this constant increased w and avoid to generate new complex root for next 15 point. Anyway it can solve offset very good but some points (nearly 15 point, but after changing constant value the point count is reduced to 5-6). I will continue to reduce unsolvable points on offsets. Is there any way to eliminate them in real-time. I think I can eliminate them by using interpolation methods, Do you think can interpolation or extrapolation methods help me to reduce also in real-time?

Below shows that how complex roots can be calculated in error(t) function. Complex roots in error(t) always form as $a+ib$ and $a-ib$. $C_1$ and $C_2$ can be calculated with using $W(t)$ and $W'(t)$. I know this way is for our differential equation but multipling with a constant is eliminate many places has problem.

$$W(t) = e^{at} [C_1 \cos(bt)+ C_2 \sin(bt)]$$

$$W'(t) = a e^{at} [-bC_1\sin(bt)+ bC_2\cos(bt)]$$

I know $W(t_0)$ and $W'(t_0)$ in the above equation.

method 2 http://ycengiz.files.wordpress.com/2013/02/method1-error.png

Anyway if you think this way is logic, I will give more details about error(t) function.

I am also adding the omega data(red dotted) with g data(green dotted) and desired g data(blue), to understand how omega and g data are affected when multiplied omega from error(t) or previous omega (it gives almost same result) with a constant value which reduces omega when complex roots has occurred. here constant is 0.94.

method 2 http://ycengiz.files.wordpress.com/2013/02/method1-error2.png

Let me know if you need more data than provided here. I haven't look your observations yet on the DE but I intend to compare soon. Thanks again.

  • $\begingroup$ I formatted the formulas in your post to improve readability. Please note that on SE answers should really be answers, not additional information about the question. Such information should be added into the question body by clicking edit under the question. One of several reasons for this requirement is that the answers normally appear in order of votes, not chronologically. This makes the discussion rather hard to follow. // I suggest that you move the content of this post into the body of the question, separating the parts with formatting such as --- (horizontal line). $\endgroup$ – user53153 Feb 26 '13 at 16:23
  • $\begingroup$ Thank you for attention and help. You are right but I think this can be a answer of my question. It includes not only more detail about question also suggests a different solution. Thank you for advises again. $\endgroup$ – yakupc Feb 26 '13 at 16:52
  • $\begingroup$ So, it looks like we are talking about the behavior on quadratic splines or cosine waves here on the way down, right? If you can give me an idea of what the ranges of the parameters in your parabolic splines can possibly be, we can figure out what is the theoretical limit for the precisions on the offsets in that range and how to (almost) match it. I'll comment more on it later when I come home. :) Thanks for the pictures. I'll look at your algorithm too but, again, later in the day. $\endgroup$ – fedja Feb 26 '13 at 17:47
  • 1
    $\begingroup$ OK, I spent a couple of days on this. I'll try to make a small writeup and post it either here or on AoPS (to make the pictures there is a piece of cake and here even posting code is a headache). The conclusion is that one can design a family of controls $K(A,V)$ applicable to the restricted motions $\dot H^2/H^{5/2}\le V,|\ddot H|/H^{3/2}\le A$ that are nearly optimal in that range (say, 30% worse than the best theoretical precision). I'm not sure what you are doing when killing $dt$. Looks like you just ignore the tangential acceleration this way, which defeats the whole idea... $\endgroup$ – fedja Mar 2 '13 at 6:38
  • 1
    $\begingroup$ OK, I posted a few pictures: artofproblemsolving.com/Forum/viewtopic.php?f=224&t=523133 You can also get the code by clicking on the images. The first one is the table of theoretical limits, the second one is the $K(2,1.6)$ controller in its normal range, and the last one is a test curve of your type (pretty harsh, actually, for this controller). The magenta line is the "trivial controller" $H^{1/4}$ and the green line is $K(2,1.6)$. I amplified the deviations 5 times to show them better using the same colors. The end of the curve is a crash test: nothing can survive... $\endgroup$ – fedja Mar 2 '13 at 23:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.