Differential equation as a function of dependent variable in the equation

Question

$$ \left\{ \begin{array}{lll} \frac{d \beta}{d t} &=& \left\{ A \beta + B(\cos{\beta}-1) \right\}^{1/2} \\ \beta(t=0) &=& 0 \end{array} \right. $$

• $A = \frac{6l}{ML^2} (F_g +F_i) \approx 3 \times 10^5$
• $B=\frac{6kl^2}{Ml^2}\approx 3.3 \times 10^5$

Here 't' stands for time and 'beta' for opening angle of structure similar to cantilever beam. A & B are constants, approximate values of which are mentioned in the image, but you can take them as any simple value for solving.

I was trying to solve a differential equation(in the image), but the problem is I cannot solve it by usual integration method and when I solve it using the runge kutta method, I get all zeros, that's because of the boundary condition, which I can judge. But when I change the boundary condition to $$ \beta(t=0) = 0.0174\frac{\pi}{180} $$ I get a linear increasing curve but actually what I must be getting is a fluctuating curve something similar to the sine wave, can you please spot any mathematical error and guide me through this? Thank you.enter image description here

Answer 1

For small values of $β$ your equation is approximately $$ \dot β=\sqrt{Aβ} $$ which is not uniquely solvable for $β(0)=0$. Obviously, this equation as well as the original equation have the zero solution. But there is also the largest possible solution of the approximation $$ β(t)=\frac{A}4t^2. $$ You can also get any solution in-between as $β(t)=\frac{A}4\max(0,t-c)^2$.

So to get a non-zero solution you can choose some $\tau$ with $Aτ\approx 0$ and start with the initial condition $β(τ)=\frac{A}4τ^2$. This solution will be monotonously increasing until hitting the first positive root of the right side, if there is one, from whereon it will stay constant at that root value.

Answer 2

You have at least $2$ problems here. The first is that you have reduced the problem to quadrature and can't handle the quadrature. The differential equation $$\frac{d\beta}{dt}=\sqrt{A\beta+B(\cos\beta-1)}$$ Could be rewritten as a quadrature $$\int_{\beta_0}^{\beta}\frac{du}{\sqrt{Au+B(\cos u-1)}}=\int_{t_0}^tdv$$ With $\beta_0=t_0=0$ this allows you to compute $t(\beta)$. If you need $\beta(t)$ you can do a couple of rounds of Newton-Raphson iterations to hit the value of $\beta$ that yields the desired $t$. Some methods of numerical quadrature like Simpson's Rule or Romberg's Rule or Clenshaw-Curtis quadrature require the value of the function at the endpoints and so aren't applicable here, but Gaussian quadrature doesn't so it should work. An even cooler method of quadrature for this problem is Gauss-Jacobi-Left Radau quadrature: it evaluates $$\int_{-1}^1(1-x)^{\alpha}(1+x)^{\beta}f(x)dx=w_0f(-1)+\sum_{i=1}^nw_if(x_i)$$ Here we would take $\alpha=0$ and $\beta=-1/2$ and let $u=\frac12\beta(x+1)$ to get $$\begin{align}t(\beta)&=\int_{-1}^1\frac{\frac12\beta\,dx}{\sqrt{\frac12\beta}\sqrt{1+x}\sqrt{A-\frac{2B\sin^2\left(\frac14\beta(x+1)\right)}{\frac12\beta(x+1)}}}\\ &=\sqrt{\frac12\beta}\int_{-1}^1(1+x)^{-1/2}\left(A-\frac{B\sin^2\left(\frac14\beta(x+1)\right)}{\frac14\beta(x+1)}\right)^{-1/2}dx\end{align}$$ And the nice thing here is that our effective $f(x)$ is well-behaved for all $x\in\mathbb{R}$ because we have factored out the singularity at $x=-1$ and we easily know $$\lim_{x\rightarrow-1^+}f(x)=\frac1{\sqrt A}$$ We need to know the nodes $x_i$ and the weights $w_0$ and $w_i$. We have computed $$w_0=2^{3/2}\left(\frac{2^{2n+1}n!(n+1)!}{(2n+2)!}\right)^2$$ For $\alpha=0$ and $\beta=-1/2$ and the nodes $x_i$ are the same as for Gauss-Jacobi quadrature for $\alpha$ and $\beta+1$. The weights $w_i$ are the same except you have to divide the Gauss-Jacobi weights $w_i$ by $(1+x_i)$ to get the Gauss-Jacobi-Left Radau weights. You can get these Gauss-Jacobi nodes and weights here: just enter the value of $n$ you want and enter $0$ for $\alpha$ and $1/2$ for $\beta$ and it computes a table which you can copy into a data file and read into your program.
Once you have the first few data points you could use any old quadrature formula or even go back to the differential equation view using Runge-Kutta because you won't be near the singularity any more.

The second problem is much more serious. Because it's a bit of a tricky point, let's look at a similar problem that we can solve. Consider the second-order differential equation of the simple harmonic oscillator: $$\ddot x=-\omega^2x$$ Let $v=\dot x$ then $$\ddot x=\frac d{dt}\frac{dx}{dt}=\frac{dv}{dt}=\frac{dv}{dx}\frac{dx}{dt}=v\frac{dv}{dx}$$ So $$v\,dv=-\omega^2x\,dx$$ And we can integrate to get $$\frac12v^2=-\frac12\omega^2x^2+C_1$$ Let the turning point be $x=A$ so that when $x=A$, $v=0$ and so $C_1=\frac12\omega^2A^2$. Solving for $v$, $$\frac{dx}{dt}=v=\pm\omega\sqrt{A^2-x^2}$$ Then start at the left turning point. We can't go backward from here so $v\ge0$ in the near future and the differential equation reads $$\frac{dx}{dt}=\omega\sqrt{A^2-x^2}$$ With the initial condition $x=-A$ when $t=0$. This has the same issues as your differential equation does, but we know how to solve it: $$\int_{-A}^x\frac{dx}{\sqrt{A^2-x^2}}=\sin^{-1}\frac xA-\sin^{-1}(-1)=\sin^{-1}\frac xA+\frac{\pi}2=\omega\int_0^tdt=\omega t$$ Then $$x=A\sin\left(\omega t-\frac{\pi}2\right)=-\cos\omega t$$ But how do we get oscillation when our first-order differential equation says that $x$ must never decrease? Well, when we reach the right turning point, $\omega t=\pi$, $x=A$, from the first-order differential equation we see that $\dot x=0$ and from the second-order equation $\ddot x=-\omega^2A<0$ so in the near future we must take the minus sign and $$\frac{dx}{dt}=-\omega\sqrt{A^2-x^2}$$ and working it through we see that the solution $x=-a\cos\omega t$ remains valid for all $t\ge0$.

Let's apply the knowledge we obtained from the simple harmonic oscillator to the problem at hand. We have $$\frac{d\beta}{dt}=\sqrt{A\beta+B(\cos\beta-1)}$$ Since $\beta=0$ when $t=0$ we are again starting from a left turning point. However if you plot $A\beta$ and $B(1-\cos\beta)$ vs. $\beta$ you will see that $A\beta>B(1-\cos\beta)$ for $\beta>0$ unless $A/B<0.7246113537767090$ or so but for this problem $A/B\approx0.9091$ so there is no right turning point and so no possibility of oscillation. What went wrong? Well, our original second-order differential equation must have read $$\frac{d^2\beta}{dt^2}=\frac12A-\frac12B\sin\beta$$ Since $0<A<b$ let $A=B\sin\beta_0$ where $\beta=\beta_0$ is the equilibrium point and then let $\beta=\beta_0+\gamma$. Then $$\frac{d^2\gamma}{dt^2}=\frac12B\left(\sin\beta_0-\sin(\beta_0+\gamma)\right)$$ So that the first integral is $$\frac12\left(\frac{d\gamma}{dt}\right)^2=\frac12B\left(\gamma\sin\beta_0+\cos(\beta_0+\gamma)\right)+C_1$$ Let the left turning point be at $\gamma=-\Gamma$ so in the near future $$\frac{d\gamma}{dt}=\sqrt B\sqrt{(\gamma+\Gamma)\sin\beta_0+\cos(\beta_0+\gamma)-\cos(\beta_0-\Gamma)}$$ If $\Gamma=\pi+2\beta_0$ then $\gamma=-\Gamma$ would be a double root of the radicand so $$\lim_{t\rightarrow-\infty}\gamma(t)=-\Gamma$$ Which isn't oscillation. Similarly $\gamma=\pi-2\beta_0$ would be a double root if $$(\pi-2\beta_0+\Gamma)\sin\beta_0-\cos\beta_0-\cos(\beta_0-\Gamma)=0$$ Which, for $A=3.0\times10^5$, $B=3.3\times10^5$, $\beta_0=1.1410966606434700$ means $\Gamma=0.4379521952834380$ when $$\lim_{t\rightarrow\infty}\gamma(t)=\pi-2\beta_0$$ Which again isn't oscillatory. The second result puts an upper bound on $\Gamma$ for oscillatory solutions. Comparing with your first integral you have $$(\beta_0-\Gamma)\sin\beta_0+\cos(\beta_0-\Gamma)=1$$ Which corresponds to $\Gamma=\beta_0=1.1410966606434700$ which is way out of range for oscillatory solutions.

All that having been said, I think it would be much easier to directly solve the original second-order differential equation numerically. At least the above analysis may give you some ideas about your initial values.

EDIT: Just for fun I made some phase space plot of the system for a range of turning points $-\Gamma$. Matlab code: Oscillator.m

% Oscillator.m

clear all;
close all;

global beta0 B

A = 3.0e5;
B = 3.3e5;
beta0 = asin(A/B);
G = 0.4;
for k = 1:10,
    y = (pi-2*beta0+G)*sin(beta0)-cos(beta0)-cos(beta0-G);
    yp = sin(beta0)-sin(beta0-G);
    err = y/yp;
    G = G-err;
    if abs(err) < 1.0e-10,
        break;
    end
end
Gmax = G;
fprintf('Gmax = %.16f\n',Gmax);
npts = 100;
legends = [];
figure;
for k = 1:npts+1,
    G = Gmax*k/(npts+2);
    options = odeset('Events',@events1,'MaxStep',2*pi/sqrt(B*cos(beta0))/100);
    [t,g] = ode45(@f,[0 0.04],[-G 0],options);
    tau(k) = t(end);
    if mod(k,10) == 1,
        plot(beta0+g(:,1),g(:,2));
        legends = [legends {['\Gamma = ' num2str(G)]}];
        hold on;
    end
end
title('Phase Space Plots for System');
xlabel('\beta');
ylabel('$\frac{d\beta}{dt}$','Interpreter','latex');
legend(legends,'Location','eastoutside');
figure;
plot([1:npts+1]*Gmax/(npts+2),tau);;
title('Period of System');
xlabel('\Gamma');
ylabel('T');

f.m

% f.m

function yprime = f(t,y);

global beta0 B

yprime = [y(2);
    B*(sin(beta0)-sin(y(1)+beta0))/2];

events1.m

% events1.m

function [value,isterminal,direction] = events1(t,y);

value = y(2);
isterminal = 1;
direction = 1;

Phase space plot:

Period as function of the parameter $\Gamma$:
As can be seen the period $T$ approaches the small amplitude limit of $$T=\frac{2\pi}{\sqrt{\frac12B\cos\beta_0}}=0.023965109359898$$ For small $\Gamma$ and increases without bound as $\Gamma\rightarrow\Gamma_{\max}=0.4379521952834382$