Numerically Solving a Second Order Differential Equation Using C

[Menninger]

Homework Statement

The Second Order Differential Equation is:
x''-u(b^2 + x^2)x'+x=0

Initial Conditions are:
x(0)=1
x'(0)=0

It is to be numerically solved for 0<=t<=500. The specific numerical method to be used isn't specified, but it must be programmed into c.

As a means to check the answer, the problem specifies that for b=1 and u<<b, the solution is approximately: (1) x(t)=2cos(t), y(t)=-2sin(t).

The Attempt at a Solution

I broke the second order equation up into two first order equations:
x'=x'
x''=u(b^2 + x^2)x'-x

Then I attempted to use Euler's Method to solve the system using the following algorithm:

(2)
x'(t+Δt)=x'+(u(b^2 + x^2)x'-x)*(Δt)
x(t+Δt)=x+x'*Δt

using the following code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char** argv) {
double t, s, y, i; /* t is the time up to which the equation is to be approximated for
* s is the number of steps to be used to calculate the approximation
* y is time length of each step
* i is time concurrent with the beginning of each step */
double x, x1, x1a, b, u, x4, x5; /* x is the value of x at the beginning each step
* x1 is the value the derivate of x at each step
* x1a is a placeholder that holds onto to the value of x1
* at the end of each step until the initial value has been
* used in the calculation of the x value at the end of the step */
x=1;
x1=0;
printf("Enter the initial value of b: \n");
scanf("%lf",&b); /* Here the user defines the initial value of b */
printf("Enter the initial value of u: \n");
scanf("%lf",&u);/* Here the user defines the initial value of u */
printf("Enter the time up to which you wish to approximate the solution of the differential equation \n");
scanf("%lf",&t); /* Here the user defines the time up to which a solution is to be approximated */
printf("%f \n", t);
printf("Enter the number of steps you wish to use to produce the approximation: \n");
scanf("%lf",&s); /* Here the user defines the number of steps which are be used to produce the approximation*/
printf("%f \n", s);
y=t/s; /* Here the user defines the temporal duration of each step */
printf("%f \n", y);
printf("First Derivative: %f Value of Function: %f\n", x1, x);
for (i=0; i<=t; i=i+y) { /* This calculates the time at the beginning of each step,
* and it terminates the loop when the final time is reached */
printf("Time: %f First Derivative: %f %f Value of Function: %f %f\n", i, x1, x4, x, x5); /* This provides the approximated position x and its derivative at the beginning and end of each step */
x5=2*cos(i);
x4=-2*sin(i);
x1a=x1+(u*(pow(b, 2)-pow(x, 2))*x1-x)*y; /* This approximates the value of the derivative of x at the beginning of the next step */
x=x+x1*y; /* This approximates the a value of x at the beginning of the next step */
x1=x1a; /* This initializes the new derivate of x for the next step */
}
return(EXIT_SUCCESS);
}

In my program (1) & (2) match in terms of periodicity, but (2) spirals outward.

Am I using the wrong technique? Do I need something more precise? And does my code accurately implement my algorithm?

 

[Mark44]

When you post code, it's considered good manners to put a [ code ] at the beginning and a [ /code ] tag at the end (without the extra spaces). I have done that, below, to make your code more readable.

Homework Statement

The Second Order Differential Equation is:
x''-u(b^2 + x^2)x'+x=0

Initial Conditions are:
x(0)=1
x'(0)=0

It is to be numerically solved for 0<=t<=500. The specific numerical method to be used isn't specified, but it must be programmed into c.

As a means to check the answer, the problem specifies that for b=1 and u<<b, the solution is approximately: (1) x(t)=2cos(t), y(t)=-2sin(t).

The Attempt at a Solution

I broke the second order equation up into two first order equations:
x'=x'
x''=u(b^2 + x^2)x'-x

The second equation above is not a first order DE. Also, your first equation is trivially true.

Your system of two first order equations should be the following:

y = x'
y' = u(b^2 + x^2)x' - x = u(b^2 + x^2)y - x

Then I attempted to use Euler's Method to solve the system using the following algorithm:

(2)
x'(t+Δt)=x'+(u(b^2 + x^2)x'-x)*(Δt)
x(t+Δt)=x+x'*Δt

using the following code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char** argv) {
    double t, s, y, i; /* t is the time up to which the equation is to be approximated for
                        * s is the number of steps to be used to calculate the approximation
                        * y is time length of each step
                        * i is time concurrent with the beginning of each step  */
    double x, x1, x1a, b, u, x4, x5; /* x is the value of  x at the beginning each step
                        * x1 is the value the derivate of x at each step
                        * x1a is a placeholder that holds onto to the value of x1
                        * at the end of each step until the initial value has been 
                        * used in the calculation of the x value at the end of the step */ 
    x=1;
    x1=0;
    printf("Enter the initial value of b: \n");
    scanf("%lf",&b); /* Here the user defines the initial value of b */
    printf("Enter the initial value of u: \n");
    scanf("%lf",&u);/* Here the user defines the initial value of u */
    printf("Enter the time up to which you wish to approximate the solution of the differential equation \n");
    scanf("%lf",&t); /* Here the user defines the time up to which a solution is to be approximated */
    printf("%f \n", t);
    printf("Enter the number of steps you wish to use to produce the approximation: \n");
    scanf("%lf",&s); /* Here the user defines the number of steps which are be used to produce the approximation*/
    printf("%f \n", s);
    y=t/s; /* Here the user defines the temporal duration of each step */ 
    printf("%f \n", y);
    printf("First Derivative: %f Value of Function: %f\n", x1, x);
    for (i=0; i<=t; i=i+y) { /* This calculates the time at the beginning of each step,
                              * and it terminates the loop when the final time is reached */ 
        printf("Time: %f First Derivative: %f %f Value of Function: %f %f\n", i, x1, x4, x, x5); /* This provides the approximated position x and its derivative at the beginning and end of each step */
        x5=2*cos(i);
        x4=-2*sin(i);
        x1a=x1+(u*(pow(b, 2)-pow(x, 2))*x1-x)*y; /* This approximates the value of the derivative of x at the beginning of the next step */ 
        x=x+x1*y; /* This approximates the a value of x at the beginning of the next step */
        x1=x1a;   /* This initializes the new derivate of x for the next step */
    }
    return(EXIT_SUCCESS);
}

In my program (1) & (2) match in terms of periodicity, but (2) spirals outward.

Am I using the wrong technique? Do I need something more precise? And does my code accurately implement my algorithm?

I didn't check that your implementation reflects the system - it wasn't clear to me that you were actually solving the right system (see my comments above).

Euler's method isn't all that precise, so even if your implementation is correct, it could that this method isn't producing accurate results once you get very far from your initial conditions. A method that is probably more precise is Runge-Kutta.

Programming tip: When you're squaring a number, it's more efficient to use multiplication than to use pow(). You won't notice a difference in your program, but it's something to keep in mind.

Word tip: there is no such word in English as "derivate".

 

[Menninger]

I only made this mistake when presenting the problem here, and not in the code itself, but when presenting the differential equation here, I miswrote it as: x''-u(b^2 + x^2)x'+x=0 .

Actually it is: x''-u(b^2 - x^2)x'+x=0

For additional clarity, the problem is presented as follows:

However, my output is a circle about half that radius, that slowly spirals outwards. Setting u to 0 produces an indefinite circle of half the desired radius that otherwise is identical to the given steady state solution. Maybe the teacher miswrote the steady state solution?

When you post code, it's considered good manners to put a [ code ] at the beginning and a [ /code ] tag at the end (without the extra spaces). I have done that, below, to make your code more readable.

Will do, thanks!

The second equation above is not a first order DE. Also, your first equation is trivially true.

Your system of two first order equations should be the following:

y = x'
y' = u(b^2 + x^2)x' - x = u(b^2 + x^2)y - x

Thanks! However, best I can tell this doesn't alter my algorithm, it just changes the variable names.

I changed my algorithm to the following to attempt to see if it made a difference:

for (i=0; i<=t; i=i+h) { 
/* This calculates the time at the beginning of each step,
                              * and it terminates the loop when the final time is reached */ 
        printf("Time: %f First Derivative: %f %f Value of Function: %f %f\n", i, x1, x5, x, x4); 
/* This provides the approximated position x and its derivative at the beginning and end of each step */ 
        x4=2*cos(i);
        x5=-2*sin(i);
        y=x1;
        y1=(u*(b*b-x*x)*y-x); 
        y=y+y1*h;
        x=x+y*h;
        x1=y;
    }

However the output was the same.

Euler's method isn't all that precise, so even if your implementation is correct, it could that this method isn't producing accurate results once you get very far from your initial conditions. A method that is probably more precise is Runge-Kutta.

Thanks! Can you link me to what the algorithm should be in the event of there being changing variables other than y and t when using Runge-Kutta? I can do 4th Order to solve a single First Order DEQ. I am not sure how to do it when dealing with a system.

At anyrate, I tried using black box Runge-Kutta software, and the output was the same as my code produced.

The radius of the circle is equal to the distance of the starting point from the origin (I think this is the right answer for part b, along with that a positive u will always cause the solution to spiral away from the steady state analog). So, as shown above, my teacher's solution works if you change the initial condition to a distance of 2.

Programming tip: When you're squaring a number, it's more efficient to use multiplication than to use pow(). You won't notice a difference in your program, but it's something to keep in mind.

Word tip: there is no such word in English as "derivate".

Thanks!

 

[Mark44]

Sorry, I haven't done anything with numerical solutions to diff. equations for a long time, so I can't offer any guidance.