Is my Diff. Eq. plot correct? d2x/dt2.

[SophiaCristina]

So, I'm relatively new on c++ and maths, don't laugh at me :(, but i made a derivative function and a d^2/dt^2 Diff. Eq. function for the derivative function, it works like this:

double Derivative(double x, double Freq)
{
 double Limit = 0.000000001;
 return ((MiniForm(x + Limit, Freq) - MiniForm(x, Freq)) / Limit); // MiniForm = sin(x * Freq) in this example
}

double DiffEq(double x, double Freq)
{
	double Limit = 0.000001;
	double Eq = (Derivative(x + Limit, Freq) - Derivative(x, Freq)) / Limit;
	return(Eq);
}

The plot i got from this is strange, this is the picture, from 0 to 3 * pi, sin(x) is the MiniForm(x, Freq), cos(x) is the Derivative(x, Freq), and this one that looks sprayed is DiffEq(x, Freq):
https://i.imgur.com/crP3vAR.png

This don't looks right, don't? Is that a proper function?

[ne555]

you have one source of error in the approximation method, using just one point in a forward derivative is O(Limit)
so the error decrease with Limit

but another source of error is in the operations:
- the sum x+Limit, where x may be as big as 10, and Limit as small as 1e-9
- the subtraction f(x+Limit) - f(x), two numbers near identical
to fix this try to increase the order of Limit, for example, Limit = 1e-6


about your graph, ¿what's the purpose of the colour? doesn't seem to have any meaning and it's quite distracting.
also, you may plot the error in your approximation

Derivative(x, 1) - cos(x)
DiffEq(x, 1) + sin(x)

[lastchance]

I tried your routines and the problem is (sort of) along the lines indicated by @ne555 - you are subtracting two numbers which are near identical. The problem is most serious in the second derivative because you are then using a discrete approximation on top of a discrete approximation.

If, in the second derivative, you change your line
double Eq = (Derivative(x + Limit, Freq) - Derivative(x, Freq)) / Limit;
to the more common finite-difference approximation for second derivatives (which doesn't rely on a double level of subtracting near-equal numbers):
double Eq = (MiniForm(x+Limit,Freq)-2*MiniForm(x,Freq)+MiniForm(x-Limit,Freq)) / (Limit*Limit );
then you will get a better result. Note that it is using the original function, not relying on the already-approximated first derivative. See code below and: https://imgur.com/fRUlTh6



As @ne55 also pointed out, you will also get a better order of accuracy if you use a symmetric difference in your first derivative; i.e. replace
return ((MiniForm(x + Limit, Freq) - MiniForm(x, Freq)) / Limit);
by
return ((MiniForm(x + 0.5 * Limit, Freq) - MiniForm(x - 0.5 * Limit, Freq)) / Limit);
There are good physical reasons for using one-sided differences in CFD, but not for what you are doing here.


Incidentally, it is a tad irritating to have to pass your natural frequency Freq through the derivative routines to get it to your function MiniForm.

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;


double MiniForm( double x, double omega )         // Guess
{
   return sin( omega * x );
}


//====== YOUR ROUTINES (AMENDED) ========

double Derivative(double x, double Freq)
{
   double Limit = 1.0e-9;
   // return ((MiniForm(x + Limit, Freq) - MiniForm(x, Freq)) / Limit);                    // Original
   return ((MiniForm(x + 0.5 * Limit, Freq) - MiniForm(x - 0.5 * Limit, Freq)) / Limit);   // Symmetric version
}


double DiffEq(double x, double Freq)
{
        double Limit = 1.0e-6;
//      double Eq = (Derivative(x + Limit, Freq) - Derivative(x, Freq)) / Limit;        // Original
        double Eq = (MiniForm(x + Limit, Freq) - 2 * MiniForm(x, Freq) + MiniForm(x - Limit, Freq)) / (Limit * Limit );   // Modified
        return(Eq);
}

//====== END OF YOUR ROUTINES ========


int main()
{
   const double PI = 3.14159265359;
   int ncycles = 2, ndx = 128;
   double omega = 1.0;
   double dx = ( 2.0 * PI / omega ) / ndx;

   ofstream out( "output.txt" );
   for ( int i = 0; i <= ncycles * ndx; i++ )
   {
      double x = i * dx;
      out << x << "  " << MiniForm( x, omega ) << "  " << Derivative( x, omega ) << "  " << DiffEq( x, omega ) << '\n';
   }
}

Plotted in gnuplot with:

p "output.txt" u 1:2 w l title "f(x)", "output.txt" u 1:3 title "dfdx", "output.txt" u 1:4 w l title "d^2 f/dx^2"

[MikeStgt]

To get the first derivative I suggest not to use only two points. How about
f'(x) = (-11f(x) + 18f(x+h) - 9f(x+2h) + 2f(x+3h))/6h where h is the current step size. Iterate it decreasing h until the result indicates "h too small".
This is what I used to use on a pocket calculator (PPCROM with HP41). Look for Numerical Differentiation.

[SophiaCristina]

Thanks everybody, @ne555, @lastchance and @MikeStgt!
@ne555 thanks, @lastchance gave me the code based on this limit that you also recommended, i used and it was correct, i also noticed that if i had used 0.0001 it would work, but the methods of you guys is more adequate and accurate...
Also, the color is bad for sin(x), but when i plot a polar function until a big number * pi, it makes a lot of knots, the beginning and the end is red so i know where is the mid point (cyan) or not, ex.: https://i.imgur.com/LKBUckh.png...

@latchance, again you helped me with the code itself, it is working perfectly (https://i.imgur.com/soX4c0q.png), thank you!

@MikeStgt, thanks for the suggestion, i never saw a derivative formula like that, i tested in wolfram alpha and it gave cos(x), but I'll use @lastchance code, it worked fine, but i'll take note of that formula, i'll be honest that i don't understood well those sine multiplications, but i'll study it.

Again, thanks everybody!

[MikeStgt]

i'll be honest that i don't understood well those sine multiplications, but i'll study it.

If possible, use calculus to create a program to compute values of the derivative. If this is not possible, there exist several polynomial approximations (that stronly smell like Taylor series), see Handbook of Mathematical Functions, table 25.2, p. 914 here -- http://www.jonsson.eu/resources/hmf/aas_hmf_master.pdf and ask a reliable math prof near you (or pick one of the lecturers you find in the Internet talking about Numerical Differentiation).

[SophiaCristina]

@MikeStgt
I'll check it out if i can, i'm not much good at maths yet, yet i can understand the concepts, but the function is working finely at least, also that is a very interesting book, thank you!
I saw on wikipedia, it have a nice article about that, i'll study it out: https://en.wikipedia.org/wiki/Numerical_differentiation

Also i found @lastchance formula here:
https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences

[MikeStgt]

Thank you for the links. There are referenced the coefficients of the formula
f'(x) = (-11f(x) + 18f(x+h) - 9f(x+2h) + 2f(x+3h))/6h, see table at https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_finite_difference , Derivative 1, Accuracy 3.
The "forward finite" method uses only f(x+nh) with n=0, 1, 2,... but no f(x-nh) what yields results even near a singularity or at discontinuity. With +h or -h you choose from which side you approach, here +0 and -0 may make a big difference. Simple example: y = abs(x), f'(+0) = 1, f'(-0) = -1.

One more word about Iterate it decreasing h until the result indicates "h too small". For successively smaller |h| values, the sequence of f' estimates should be monotonic and |f'_i+1 - f'_i| should be monotonic and decreasing. If either condition is violated |h| is probably so small that it causes errors.

But that is not of much help in your a. m. case.

[SophiaCristina]

Hmm, now i can see where you got those multiplications, you got this table function and multiplied each fraction by 6... Interesting, numerical analysis is something i should give more focus upon it...
Interesting example too, but, why "-0" gives "-1" if it is "abs(x)"?

It is helpful in matter of information and if someone else check this page, also that is basically what happened, h was too small between both derivatives, instead of using a formula to another formula, i should have used a single second-order formula...
Don't worry, i ignore nothing, i'll take those tips in mind...

[MikeStgt]

but, why "-0" gives "-1" if it is "abs(x)"?

Not abs(x) but its first derivative, abs(x)d/dx at x=0-0.
In case you had to program the abs function it could look like

if x > 0 return x;  // y = x (y is positive as x)
else return -x;     // y = -x (make y positive) 

If you picture the graph it consists of two diagonals, for x<=0 (below 0 or "to the left of x=0") it is a diagonal line in the 2nd quadrant from the upper left to the origin (0, 0). For x>0 this is mirrored at the ordinate axis. The lines' slope (first derivative) is -1 for all x<0 and 1 for all x>0. And what at precisely x=0? Numerically no chance, you have to compute at x=0 (abs(x+h)-abs(x))/h or from the lhs (abs(x-h)-abs(x))/-h. In the first case you get h/h=1, in the second case h/-h=-1 (BTW, using a mid-point formula the result would be 0 what is wrong).
In calculus the direction of approach to a function's singularity or discontinuity at x is indicated by x+0 or x-0. That's all. But enough for a smile for mathematicians when they can show their students that a 0, null, zero, nada with a plus or minus sign may make a big difference.

[lastchance]

The mathematical function abs(x) is simply NOT differentiable at x = 0. It is nonsense to suggest that it has "different derivatives from each side". A function has at most ONE value at a point - by definition. EVERY difference formula for dydx would get it wrong for y = abs(x) ... because that function simply doesn't have a derivative at x=0.

Central-differenced formulae are inherently more efficient than one-sided differences. A one-sided difference formula making 4 function evaluations would have order 3, whereas a centrally-differenced formula making 4 function evaluations would have order 4. For a given number of function evaluations a centrally-differenced formula is more accurate.

There are some reasons for using one-sided differences for advective derivatives in CFD, because information travels in the direction of fluid flow - and you would use a backward difference in that case, not a forward one. Just about everywhere else a centrally-differenced formula would be far preferable.

Apart from anything else, the OP's original question concerned the second derivative, not the first. The rest of this thread has been highly irrelevant.

[MikeStgt]

The mathematical function abs(x) is simply NOT differentiable at x = 0.

D'accord. Agreed. (Are you really sure?)
Now -- little subject drift -- how do you determine if (typical example) sin(x)/x has a removable singularity at x=0? I used abs(x) only as a not so nice example to explain x+0 and x-0.

A function has at most ONE value at a point - by definition.

By definition sounds like by agreement. What about sqrt(x)?

Edit: see strike through

[Ganado]

What about sqrt? What are you implying? sqrt(x) is a function both mathematically and by C++.

sin(x)/x has a removable singularity because its limit from the left matches its limit from the right (and both exist).

[lastchance]

sqrt(x) isn't differentiable at x=0 either.

It is a function with two branches: if you want it to be continuous, let alone differentiable, you will have to stay on one of those branches (as c++ does).

@MikeStgt, I have no idea why you are straying into removable singularities. Once again, this has nothing to do with the OP's thread, or finite differences. Any mathematically "pathological" point could potentially cause any particular numerical approximation to fail, and shouldn't be used as a diversion. "Subject drift" of this size would be better dealt with by a separate thread.

[MikeStgt]

I have no idea why you are straying into removable singularities

I assumed you could know about the prerequisites for singularities to be removable. The subject drift started before, when I tried to explain why I suggest a 'single sided' formula even the errors of central-differenced formulae are smaller.
http://www.cplusplus.com/forum/beginner/254949/#msg1118838
(Also see last phrase there.)

BTW, two appends before that I recommended "If possible, use calculus to create a program to compute values of the derivative." Almost no errors to consider.

As the subject drift is strongly related to the course of discussions here, I see no need to open a new thread. At least I do not have any questions in this respect.

[SophiaCristina]

Just wanted to say i have not ignored the replies, but i don't know much what to reply, you guys know way more...

@MikeStgt, I can see what you mean now on abs(x):
https://www.wolframalpha.com/input/?i=abs(x)
In wolfram it outputs derivative as x/|x|
I can see now what you mean about some functions using "-h"...
I'll use those numerical analyzes formulas from that book when i run in to a problem...

But I also read the @lastchance concept...
I think I got what you mean, computer aren't much smart, so they will try to process abs(0) and return erratically, is that so?
I'll use in standard the CFD, also, i'm using almost for what you mentioned, i'll take the derivative of acoustic pressure so i need to take second-order derivative from "d^2pressure/dx^2" and "d^2pressure/dt^2"...

Also, can't i fix '0' problems with "if 0 then 1e-9", or "if -h then |h|"?

As you said "But that is not of much help in your a. m. case.", I'll take that in mind, when i run into that problem, i'll remember those tips... ;)