piecewise differential equations solution

Question

I have to solve the exact solutions of the following piece-wise differential equations. $$x''+0.4x'+x=1\ (x'>2.5)$$ $$x''+0.4x'+x=-1\ (x'<-2.5)$$ $$x''-0.4x'+x=0\ (-2.5<x'<2.5)$$ $$x(0)=0.1,\ x'(0)=0.0$$ The phase plot of the solutions is expected to be like the image below.

Euler Methods

The above phase plot is generated by Euler methods and is proved to be the right answer. However, when finding the exact solutions using the below C program, the result is not like the image above. With the decided initial conditions, it begins with the third equation.When time reaches 17.1s (when $x'$ reaches boundary) it jumps to the first equations but the solution $x'$ it generated immediately return it to the third equation which is not continuous like the image.

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <math.h>

#define DT 20
#define T DT*1000
#define del 0.2
#define del2 0.2
#define omega1 sqrt(1-del*del)
#define omega2 sqrt(1-del2*del2)
#define Th 2.5
void main(void)
{
FILE *fp;
int i;
\\define the initial conditions x(0)=0.1, x'(0)=y(0)=0, h=0.05
double x = 0.1, t = 0.0, h = 1.0 / DT, y = 0.0,x0,y0; 

fp = fopen("Kai.csv", "w+");
fprintf(fp, "tau, x, x'\n");


for (i = 0; i <= T; i++) {
x0 = x; \\ Save the new initial conditions for using when x' reach the boundary
y0 = y;
if (y > Th) {
    while (y > Th) {
        \\Solution for x, x' with initial conditions x0, y0
        x = exp(-del*t)*((x0 - 1)*cos(omega1*t) + (1 / omega1)*(y0 + del*(x0 - 1))*sin(omega1*t)) + 1; 
        y = exp(-del*t)*(y0*cos(omega1*t) - (1 / omega1)*(del*y0 + (x0 - 1))*sin(omega1*t));
        fprintf(fp, "%f, %f,%f\n", t, x, y);
        t = t + h;
    }
}
else if (y < -Th) {
    while (y < -Th) {
        x = exp(-del*t)*((x0 + 1)*cos(omega1*t) + (1 / omega1)*(y0 + del*(x0 + 1))*sin(omega1*t)) - 1;
        y = exp(-del*t)*(y0*cos(omega1*t) - (1 / omega1)*(del*y0 + (x0 + 1))*sin(omega1*t));
        fprintf(fp, "%f, %f,%f\n", t, x, y);
        t = t + h;
    }
}
else {
    while ((y>=-Th) && (y<=Th)) {
        x = exp(del2*t)*(x0*cos(omega2*t) + (1 / omega2)*(y0 - del2*x0)*sin(omega2*t));
        y = exp(del2*t)*(y0*cos(omega2*t) + (1 / omega2)*(del2*y0 - x0)*sin(omega2*t));
        fprintf(fp, "%f, %f,%f\n", t, x, y);
        t = t + h;
        }
    }
}
fclose(fp);
}

I can't find what is the problem with my code. Please show me where I went wrong in this code. Any advice would be appreciated.

Answer 1

Your propagator formulas for $x''+2\delta x'+x=f$, $ω=\sqrt{1-δ^2}$, for the values at $t=t_0+Δt$ should be \begin{align} x(t)&=f+e^{-δ\,Δt}\Bigl((x_0-f)·\cos(ω\,Δt)+\frac1ω(y_0 + δ·(x_0 - f))·\sin(ω\,Δt)\Bigr) \\ y(t)&=e^{-δ\, Δt}\Bigl(y_0·\cos(ω\, Δt) - \frac1ω(δ·y_0 + (x_0 - f))·\sin(ω\, Δt)\Bigr) \end{align} where $x_0=x(t_0)$ and $y_0=y(t_0)$ are the values at the last step. This is largely correct in your code, but you are using $t$ instead of the time increment $Δt=t-t_0$ leading to the observed jumps at the first change of the cases.

Easy test: for $Δt=0$ the formulas should reproduce the values $x_0,y_0$, which they do.

After correcting the code by introducing t0 along x0,y0 and using (t-t0) in the formulas, the solution is indeed spiralling outwards towards a, not exactly circular, limit cycle.