1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerically Solving a Second Order Differential Equation Using C

  1. Jan 5, 2013 #1
    1. The problem statement, all variables and given/known data

    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).

    3. 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?
     
  2. jcsd
  3. Jan 6, 2013 #2

    Mark44

    Staff: Mentor

    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.

    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
    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".
     
  4. Jan 7, 2013 #3
    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:
    2s8mnpt.jpg

    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?

    Will do, thanks!

    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:
    Code (Text):
    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.

    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.

    316um8k.jpg

    4rywza.jpg

    2sb66x2.jpg

    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.

    Thanks!
     
    Last edited: Jan 7, 2013
  5. Jan 7, 2013 #4

    Mark44

    Staff: Mentor

    Sorry, I haven't done anything with numerical solutions to diff. equations for a long time, so I can't offer any guidance.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook