Using GSL for solution of a differntial equation

[Ahzmandius]

Hi guys,

Currently I am trying to solve a differential equation in c++.

V - vc(t) - R/C * vc'(t) - L/C * vc''(t)=0
V = DC Voltage = const
vc (t) = Voltage of the Capacitor
L = Inductivity
C = Capacity
R = Resistance

It's a series connection of R, C, L and a DC Voltage supply.
To solve this differential equation in C++ I am trying to make use of the GSL ODE-Solver.
I tried to implement the example in GSL's documentation and it worked fine.
The example solves the second-order nonlinear Van der Pol oscillator equation:

u''(t) + µu'(t)(u(t)^2 - 1) + u(t) = 0

https://www.gnu.org/software/gsl/doc/html/ode-initval.html:

After this I tried to solve the equation below and it did not work for me.

V - vc(t) - R/C * vc'(t) - L/C * vc''(t)=0

It got frozen after the first print order.
I am not quite sure whether my jacobi matrix entries are correct or not.

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int
func(double t, const double y[], double f[],
    void* params)
{
    (void)(t); /* avoid unused parameter warning */
    double * mu = (double *)params;
    f[0] = y[1];                        
    f[1] = (mu[2]-y[0]-y[1]*mu[0])*1/mu[1];
    return GSL_SUCCESS;
}

int
jac(double t, const double y[], double* dfdy,
    double dfdt[], void* params)
{
    (void)(t); /* avoid unused parameter warning */
    double * mu = (double*)params;
    gsl_matrix_view dfdy_mat
        = gsl_matrix_view_array(dfdy, 2, 2);
    gsl_matrix* m = &dfdy_mat.matrix;
    gsl_matrix_set(m, 0, 0, 0.0);
    gsl_matrix_set(m, 0, 1, 1.0);
    gsl_matrix_set(m, 1, 0, -1/mu[1]);//mu[1] = CL
    gsl_matrix_set(m, 1, 1, -mu[0]/mu[1]);//mu[0] = CR
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;
    return GSL_SUCCESS;
}

int
main(void)
{
    double R = 10000;
    double L = 10 * 10 ^ -3;
    double C = 10 * 10 ^ -6;
    double mu[3] = {C*R,C*L,10};
    gsl_odeiv2_system sys = { func, jac, 2, &mu };

    gsl_odeiv2_driver* d =
        gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd,
            1e-1, 1e-1, 0.0);
    int i;
    double t = 0.0, t1 = 100.0;
    double y[2] = { 0.0, 0.0 };

    for (i = 0; i <= 100; i++)
    {
        double ti = i * t1 / 100.0;
        int status = gsl_odeiv2_driver_apply(d, &t, ti, y);

        if (status != GSL_SUCCESS)
        {
            printf("error, return value=%d\n", status);
            break;
        }

        printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

    gsl_odeiv2_driver_free(d);
    return 0;
}
```end cpp

[SeppJ]

Too late for me too look at that today, but I wanted you to know, that you are programming C. Not C++, as you claim, nor C++/NET where you asked your question. Big difference. All three are completely different languages (albeit some similar syntax and common ancestry). Might be important to know when you seek help.

[Ahzmandius]

I am not quite sure what you mean. I am working with Visual Studio for C++ (at least it was what I have installed) and had no such warnings what so ever.

[Jockelx]

@Ahzmandius

Hi,

in C (C++) '^' means xor, not pow!

[wob]

@Jockelx sagte in Using GSL for solution of a differntial equation:

@Ahzmandius

Hi,

in C (C++) '^' means xor, not pow!

Correct.
So double L = 10 * 10 ^ -3; is L = -103
you might be looking for: double L = 10e-3; which is L = 0.01

What @SeppJ wrote:
#include <stdio.h> is a C header. It is available in C++ but you usually don't use it unless for C compatibility reasons (and then you might want to use <cstdio> instead). Also, the void* params is not something you would do in idiomatic C++, neither is the cast (double *)params; (it's called "C style cast" and those are frowned upon in C++). Also, in int main(void) the void is left out in idiomatic C++. So your code might compile in C++ mode, but it doesn't look like C++ and certainly isn't something a C++ developer has written. It looks like C instead. That's why your question has been moved.

[*john 0]

Ok, it's friday so I fixed the most obvious C++ errors for you and the program does now run.

#include <cstdlib>
#include <iostream>
#include <iomanip>
//#include <cstdio>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int func(const double t [[maybe_unused]], const double y[], double f[], void* params) {
    double* mu = static_cast<double*>(params);

    f[0] = y[1];
    f[1] = (mu[2]-y[0]-y[1]*mu[0])*1/mu[1];

    return GSL_SUCCESS;
}

int jac (double t [[maybe_unused]], const double[], double* dfdy, double dfdt[], void* params) {
    double* mu = static_cast<double*>(params);

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 2, 2);
    gsl_matrix* m = &dfdy_mat.matrix;
    gsl_matrix_set(m, 0, 0, 0.0);
    gsl_matrix_set(m, 0, 1, 1.0);
    gsl_matrix_set(m, 1, 0, -1/mu[1]);//mu[1] = CL
    gsl_matrix_set(m, 1, 1, -mu[0]/mu[1]);//mu[0] = CR
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;
}

using std::cout, std::setprecision, std::setw;

int main () {
    constexpr double R = 10000;
    constexpr double L = 10e-3;
    constexpr double C = 10e-6;
    double mu[3] = {C*R,C*L,10};
    gsl_odeiv2_system sys = {func, jac, 2, &mu};

    std::unique_ptr<gsl_odeiv2_driver,decltype(&gsl_odeiv2_driver_free)>
      d(gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-1, 1e-1, 0.0), &gsl_odeiv2_driver_free);
    // gsl_odeiv2_driver* d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-1, 1e-1, 0.0);

    double t = 0.0;
    constexpr double t1 = 100.0;
    double y[2] = {0.0, 0.0};

    for (int i = 0; i <= 100; i++) {
        double ti = i * t1 / 100.0;
        int status = gsl_odeiv2_driver_apply(d.get(), &t, ti, y);

        if (status != GSL_SUCCESS) {
            cout << "error, return value=" << status << "\n";
            return EXIT_FAILURE;
        }

        cout << setw(11) << setprecision(5) << t << " " << setw(11) << setprecision(5) << y[0] << " " << setw(11) << setprecision(5) << y[1] << "\n";
        //printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

    // gsl_odeiv2_driver_free(d);

    return EXIT_SUCCESS;
}

Just a short update to use a more C++ style of resource allocation.

[SeppJ]

That does not really help. Now we are just programming C in C++. There exists an open-source C++-wrapper for GSL, which could be used to write "real" C++. I have no idea if that still works, though, it seems rather old. It should not be too hard to write your own wrapper, just a lot of busy work because GSL is rather big.

Though I doubt that @Ahzmandius even wants to write C++. Everything points to them wanting to speak C and just being confused.

[Ahzmandius]

Thank you all.

The values for R, C and L caused the issue. They got too big and so the calculation got ages to get finished.

Is there somewhere some sort of list that explains why the c++ standard is the way it is? 
I already knew that in c++ you do not write int main(void){...} but why?

For example, I found an explanation why you should not use c style cast, but it would be pretty tedious to search for every c++ standard individually.

[SeppJ]

@Ahzmandius sagte in Using GSL for solution of a differntial equation:

Thank you all.

The values for R, C and L caused the issue. They got too big and so the calculation got ages to get finished.

Is there somewhere some sort of list that explains why the c++ standard is the way it is? 
I already knew that in c++ you do not write int main(void){...} but why?

For example, I found an explanation why you should not use c style cast, but it would be pretty tedious to search for every c++ standard individually.

While there is a comprehensive documentation of the C++ standard (http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2017/n4713.pdf), that is not really what you want. That is for the guys writing compilers, and for settling disputes about exotic corner cases.

What you really want, is a good book about C++. That will help you learn the language. If the code above (which is C) was written by you, then you will have to basically re-learn everything, except the basic syntax. But you can probably skip the basics and start with one of the more advanced books, like "A tour of C++". If you mostly copy&pasted the code above, without really understanding all of it, then you should start with one of the introductory books, instead. There's no shame in that, but it is important that you are honest about your own skills.

You might also be interested in language references, which are much more readable than the language standard:
https://en.cppreference.com/w/
http://www.cplusplus.com/reference/
While these cannot teach you to speak the language fluently, they do explain all the words and grammar.

[Ahzmandius]

@SeppJ
Well, this is still pretty vague.

I am aware for example that the printf() is part of c and cout is part of c++.
But still the question remains, what are some if not all concrete issues that can appear by using printf() instead of cout besides that it's not the standard c++.

The thing is, I am not a c++ developer nor software developer in general. I am an engineer who utilizes c/c++ to solve some problems. I am not writing the code, so a purist c++ developer is pleased.

Not only that, but I also do not see a reason to make everything a class when c++ supports a procedural programming too (I guess it was what you meant by "Now we are just programming C in C++."). Correct me if I am wrong.

I mean, I know what the benefits are, but at this simple example it's like cracking a nut with a sledgehammer. I also can start splitting the declarations and definitions of every function by using header and cpp files, still the nut cracking in my eyes.

This example should solely show me what the gsl-lib is capable off and how it is working. The overhead of the c++ standard is not my scope of view at this point. When I do implement it in my mfc program I will use all the standards like const expressions, const variables in function declarations etc. At this point, it really does not matter in my eyes.

[SeppJ]

Phew. I don't think it is possible to explain all the difference between C and C++ in concrete technical terms. It's always going to be vague, unless you write a whole book on the topic (And that is just, what books such as "A Tour of C++" are!). So, let's be vague: C and C++ are both capable of the same technical results. But C++ comes with a lot of very powerful tools for generalization. It might be harder to write such general code than it is to achieve a particular result in C, but it is way easier to re-use that same code for a whole bunch of similar problems than it would be to re-use the C code.

The printf vs. cout is a good example for this: On a surface level, printf is more intuitive and easier to use. But the big difference is, that printf supports only the built-in types of C and will never be able to do anything more. Whereas you can extend cout to do pretty much anything. Which is heavily used by almost any C++ library. Because of this, you should strive to use the C++-abstractions everywhere and not rely on easy shortcuts in C. Your code will be able to easily interface with any C++ library.

Look at your GSL example: You are using the C interface. Was the experience easy and intuitive for you? You have provide GSL with functions with a very particular signature. These functions get their parameters via void-pointers that you have to cast, using your knowledge about the GSL callstack. Then you have to manually allocate resources for GSL via gsl_odeiv2_driver_alloc_y_new (so easy!), apply your functions to these resources, and you must not forget to free these resources in the end. In C++ you would have constructed some kind of GSL object and given it any kind of callable object. And that's it. But of course you need to be familiar with object- and callable-abstractions to do that. Therefore, you should always strive to use idomatic C++ in any context, even for simple examples. That is how you get familiar with these abstractions.

Unless you just want to do C, not C++. Which is legit. A lot of people are very happy with that and they have good reasons. C++ does lend itself to overthinking simple things, there is not denying that. Just don't try to mix the two languages. They are very different, especially on an advanced level. As a beginner it might seem that C++ is just C with classes, or that C is C++ without classes, but neither of these is true.

[*john 0]

@Ahzmandius sagte in Using GSL for solution of a differntial equation:

@SeppJ
Well, this is still pretty vague.

I am aware for example that the printf() is part of c and cout is part of c++.
But still the question remains, what are some if not all concrete issues that can appear by using printf() instead of cout besides that it's not the standard c++.

printf is not typesafe (It is part of the C++ standard). Therefore simple errors in the format string can easily result in undefined behavior and your program will have some weird side effects. That is not good. Well, iostream is not so pleasant to use. There are type safe replacements for printf like boost format. By the way boost is always a good place to search for solutions for C++ problems. Many of the experts provides code for the boost project. So the probability is good, that someone else solved your problem before.

The thing is, I am not a c++ developer nor software developer in general. I am an engineer who utilizes c/c++ to solve some problems. I am not writing the code, so a purist c++ developer is pleased.

The purpose of this programming style is not to please other C++ programmers, but to eliminate many potential bugs in your software. C++ is not the easiest programming language. I prefer often modern Fortran (at least Fortran 2003) for programs, because it is much easier to use for scientist and engineers.

Not only that, but I also do not see a reason to make everything a class when c++ supports a procedural programming too (I guess it was what you meant by "Now we are just programming C in C++."). Correct me if I am wrong.

There is no reason to use classes, if there are other good solutions.

This example should solely show me what the gsl-lib is capable off and how it is working. The overhead of the c++ standard is not my scope of view at this point. When I do implement it in my mfc program I will use all the standards like const expressions, const variables in function declarations etc. At this point, it really does not matter in my eyes.

If you don't know there are other good libraries for solving ODEs like SUNDIALS form the LLNL. This library has C, C++ and Fortran APIs.

[*john 0]

@SeppJ sagte in Using GSL for solution of a differntial equation:

That does not really help.

I know. I would use SUNDIALS, which has a C++ API.

[hustbaer]

@Ahzmandius sagte in Using GSL for solution of a differntial equation:

Is there somewhere some sort of list that explains why the c++ standard is the way it is?

I don't think there is.

I already knew that in c++ you do not write int main(void){...} but why?

In C, you can use int foo(); to declare a function without specifying what the parameter list is. When calling the function, the compiler will simply assume that the argument types exactly match the parameter types. Like this:

int foo();

int main(void) {
    return foo(1, 2LL, 3);
}

AFAIK this goes back to K&R C, which did not even allow to specify the parameter list in a function declaration.

Back to the example above: If foo has parameters that are compatible with int, long long, int, then this will work. It's not type safe though. If you'd call the function with foo(1, 2, 3) instead of foo(1, 2LL, 3), things would potentially stop working (or, depending on several factors, it could work "by accident").

C++ OTOH doesn't allow declaring functions without specifying the parameter list, you always have to specify the parameter list now. That means the syntax int foo(); wasn't used for anything anymore. So C++ is free to re-define what int foo(); means. I think re-defining it as declaring a function which has an empty parameter list (vs. an unknown one) makes sense. IMO it's more "natural" to write int foo(); for an empty parameter list than to write int foo(void); - and of course it's shorter.

I think that's the reason why. But like with all questions of why something is like it is, a definitive answer could only come from the person/people who made it that way. But I also think that - as with most such questions - the answer to why isn't all that interesting. Much more interesting is the answer to the question "what is it good for?". And the answer to that should be pretty easy: it's shorter and more natural.

But still the question remains, what are some if not all concrete issues that can appear by using printf() instead of cout besides that it's not the standard c++.

A big disadvantage of printf is that it's not type-safe. Stuff like printf("%d", functionThatReturnsLongLong()) compiles, but does not work as intended. And if the function name isn't functionThatReturnsLongLong but more like getFoo, then such mistakes can be pretty hard to spot. At some time in the last 2 decades though, compilers have began to issue warnings if the argument list doesn't match the format string. So it's less of an issue now than it was 10 or 20 years ago.

[*john 0]

@hustbaer sagte in Using GSL for solution of a differntial equation:

AFAIK this goes back to K&R C, which did not even allow to specify the parameter list in a function declaration.

That's not correct. K&R C declares the type of parameters differently.

do_something(s)
char[] s;
{
    double r;
    /*  some code */

    return r;
}

main()
{
    double do_something();
    double s;
    static char *str[] = "123456789";

    s = do_something(str)
}

If you want, you can download a copy of the 1978 edition of "The C Programming language" from the internet archive. Robert Nordier released a UNIX 7 Version which can be run on a virtual machine, which includes an ancient K&R C compiler.