Pointers and a problem to solve with the

You won't get exactly the same columns - you are doing a numerical approximation to a differential equation - the key word is APPROXIMATION. Do 100 steps with h=0.01 to see if it gets closer.

Do you really need to use pointers?
I've no idea whether your existing step is right, but second-order explicit RK amounts to

void step(double *pt, double *px, double h)
{
      *px = *px + h * 0.5 * (  fk(*pt, *px)  +  fk(*pt+h, *px+ h*fk(*pt, *px) ));
}

4th-order RK is much more common, but I'd hate to write that out with pointers.

Look at this snippet:

void step(double *pt, double *px, double h)
{
      double *px_2;
      px_2 = new double;
      double t_1, x, x_1, x_2, x_0 = 1;
      *px = *px + h * fk(*pt, *px);
      *px_2 = *px + h * 0.5 * (fk(*pt, *px) + fk(*pt, *px));
      delete px_2;
}

What is the point of your px_2 variable? You do some computations but never use the variable after the computations.

You also seem to have several variables that you don't use.

What is the purpose of this function?

As written it really just:

void step(double *pt, double *px, double h)
{
      *px = *px + h * fk(*pt, *px);
}

Is that what you intended?

Repeater (3046)

While we're looking at this snippet:

void step(double *pt, double *px, double h)
{
      double *px_2;
      px_2 = new double;
      double t_1, x, x_1, x_2, x_0 = 1;
      *px = *px + h * fk(*pt, *px);
      *px_2 = *px + h * 0.5 * (fk(*pt, *px) + fk(*pt, *px));
      delete px_2;
}

As well as what's the point of the px_2 variable, what's the point of creating it using new? It's better to put local variables on the stack, without new:

void step(double *pt, double *px, double h)
{
      double px_2;
      double t_1, x, x_1, x_2, x_0 = 1;
      *px = *px + h * fk(*pt, *px);
      px_2 = *px + h * 0.5 * (fk(*pt, *px) + fk(*pt, *px));
}

In this case there's no reason for px_2 to exist at all, but it illustrates using a local variable on the stack rather than going through the pointless (and dangerous) labour of new and delete when you don't have to.

Last edited on

Yeah, unfortunately the assignment says to use pointers :( so I have no choice, but to use them.

Do you mean that you dont know if my first step is right or the second? What can be wrong there? As for second-order RK amount I changed it to your proposed values, but still got the same values as I got with my version of 2nd order RK. When it's *pk_1 = ...it's the ones I described in my last post and even lower when it is *pk =... (1.000, 0.8, etc).

Fortunately, I don't have to write the 4th-order RK. My lecturer says I should be able to get the values (as posted in the post above) from using 2nd-order RK.

When I tried doing approximation for 100 times I still got the same first 10 values (0.9000 instead of 0.9050, 0.8100 and so on). Probably something is wrong with my declarations or usage of pointers. Or maybe the code might be smoother at some point. Anyhow thats what I think. Thank you very much for all the help!

jib as for the point of *px_2 it's mainly that changing *px to it gave me higher values (0.9 for t = 0.1) as compared to *px ( 0.8 for t = 0.1). To be honest I don't know pointers very well so just did trial and error to hopefully get the results closer to what I should. I know it's not the perfect approach:/ However now I noticed my mistake the void function only takes the value of *px so *px_2 is unusable. So when I changed *px instead of *px_2 and deleted first *px equation I got the same values (0.9, 0.8, etc.)

You say several variables that I don't use, which ones you mean? Thanks.

The code should basically have (where x_1 here is *px in my case, but I don't know how to write x_1~ in pointers as it is necessary for 2nd equation) two RK equations:

x_1~ = x_0 + hf(t_0, x_0); h is an interval and f is function
x_1= x_0 + h/2[f(t_0, x_0) + f(t_1, x_1~)]

Repeater, thanks for noticing it. I went through new/delete routine as without them I get segmentation fault: 11 when declaring *px_2.

When I'm trying to do equations to be equal to the one given by me where are x_1~ and x_1, I get the values of x (in the answer) to be equal to 1.0000 all 10 times. Here's the snippet code:

void step(double *pt, double *px, double h)
{
      double *px_2;
      px_2 = new double;
      px = new double;
      double x_0 = 1, t_1 = 0.1;
      *px_2 = *px + h * fk(*pt, *px);
      *px = *px + h * 0.5 * (fk(*pt, *px) + fk(t_1, *px_2));
      delete px_1;
      delete px;
}

Any ideas what I might be doing wrong?

The problem basically says this: Write 2nd order Runge Kutta program to integrate differential equation. The conditions given are: (d/dt)x = -x with initial condition x(0) = 1 and step h = 0.1. The answer I should get is shown at the top post. Program should be made of three functions (double fk(double t, double x), void step(double *pt, double *px, double h) and int main()). Hope now you have a better idea..

Last edited on

jib as for the point of *px_2 it's mainly that changing *px to it gave me higher values (0.9 for t = 0.1) as compared to *px ( 0.8 for t = 0.1).

Where does anything pertaining to px_2 change px in any way?

This calcualtion: *px_2 = *px + h * 0.5 * (fk(*pt, *px) + fk(*pt, *px)); doesn't modify px in any way.

You say several variables that I don't use, which ones you mean?

All of the following: double t_1, x, x_1, x_2, x_0 = 1;

Last edited on

jib thank you for your response. That's true px_2 doesnt modify px, it's useless. I need to work with *px and make it work.

Ah yes, those variables I just wrote when I first started writing and forgot to delete them. They're not needed.

While you are at it, you can update *pt during the step as well. At the moment you are setting it explicitly without using pointers in your main loop.

Found out my problem. It was falsely understood RK method. The snippet correctly should look like this:

void step(double *pt, double *px, double h)
{
      double k_1, k_2;
      k_1 = fk(*pt, *px);
      k_2 = fk(*pt + h, *px + h * k_1);
      *px = *px + h / 2 * (k_1 + k_2);
}

Unfortunately still don't know how to use pointers. We'll need to read on them more and hopefully do more practise.

Last edited on

Thanks lastchance! How could I update *pt in main? There is a declaration that updates value (t) or am I missing something? I put updating *pt in step() function. That might be neater.

Last edited on

No, update it in step(), if you aren't going to set it explicitly as you do currently (and not unreasonably).

Once you've introduced k_1 and k_2 there doesn't seem much point in using pointers. Maybe your teacher doesn't understand passing by reference.

Have you modified fx() from the code posted in your first post? If so how? If not why are you passing two parameters but only using one?

Unfortunately still don't know how to use pointers. We'll need to read on them more and hopefully do more practise.

Yes pointers can be tricky, which is why in C++ many people prefer to use pass by reference instead of pointers.

update it in step()

Thanks lastchance. Done it. How could I pass by reference in my code? Tired to think on my own, will go to sleep.

Have you modified fx() from the code posted in your first post? If so how? If not why are you passing two parameters but only using one?

jib I didnt modify fk() function. It's the same as in the first post. I'm passing now two parameters to step() function. Taught could get away with passing just one, but after reading the problem it implicitly stated that both parameters should be passed as shown in the updated code below.

#include <iostream>
#include <stdio.h>
#include <cmath>

void step(double *pt, double *px, double h);

double fk(double t, double x);

int main()
{
      double t = 0, x = 1, h = 0.1;
      printf("Equation dx/dt = -x solution when x(0) = 1 \n \n t    x        exp(-t) \n \n");
      for(int i = 0; i <= 10; i++)
      {
              printf("%3.1lf  %5.4lf  %5.4lf \n", t, x, exp(-t));
              step(&t, &x, h);
      }
}

void step(double *pt, double *px, double h)
{
      double k_1, k_2;
      *pt = *pt + h;
      k_1 = fk(*pt, *px);
      k_2 = fk(*pt + h, *px + h * k_1);
      *px = *px + h / 2 * (k_1 + k_2);
}

double fk(double t, double x)
{
      return -x;
}

Last edited on

Taught could get away with passing just one, but after reading the problem it implicitly stated that both parameters should be passed as shown in the updated code below.

Then you misread something. Either you can pass only one parameter or you need to use both parameters in that function for something. There is no reason to pass a parameter you never use into the function. In fact the way your function is written there is no real need for it anyway. Something like k_1 = -1.0 * (*pt, *px); would be all you need to do.

Lastly you really should be using meaningful variable and function names to help document the purpose of the variables and functions.

How could I pass by reference in my code?

Just use
void step( double &t, double &x, double h )
The use of the &t and &x means that you can refer to them in step as just t and x, and their changed values will be passed back to the calling routine (here, main()).

#include <stdio.h>
#include <cmath>

void step( double &t, double &x, double h );
double fk( double t, double x );

int main()
{
   double t = 0, x = 1;                     // Initial boundary condition
   double h = 0.1;                          // Step size
   printf( "Equation dx/dt = -x, boundary condition x(0) = 1\n\n" );
   printf( "t    x       exp(-t)\n\n");
   for ( int i = 0; i <= 10; i++ )
   {
      printf( "%3.1lf  %5.4lf  %5.4lf \n", t, x, exp( -t ) );
      step( t, x, h );
   }
}


void step( double &t, double &x, double h )
{
   double k1 = fk( t, x );                  // First estimate of slope (Forward Euler)
   double k2 = fk( t + h, x + h * k1 );     // Second estimate of slope

   x = x + 0.5 * h * ( k1 + k2 );           // 2nd-order explicit Runge-Kutta; aka Heun's method
   t = t + h;                               // Update time
}


double fk( double t, double x )
{                                           // code up whatever equation is to be solved
   return -x;                               //    dx/dt = f( t, x )
}

If you leave fk( double t, double x ) as it is, with two parameters, then you have a flexible program that will solve (in principle) more complicated problems; e.g.
dx/dt = t - x
with minimal change. If fk was only a function of x then it would be pointless to have a differential equation solver, because you could just separate variables:
dx/f(x) = dt
and integrate to get a solution.

I say in principle because you might like to try solving
dx/dt = -40x
instead. It has the positively benign exact solution exp(-40t), but see what happens if you keep the same step size.

BTW, line 23 in your code should be AFTER the setting of k_1 and k_2. If fk() had been a function of time t as well (no, I know it isn't here) then the time level would be wrong.

Last edited on

closed account (48T7M4Gy)

#include <iostream>
#include <cmath>

using namespace std;

void step(double *, double *, const double h);

double f( const double );

int main()
{
    double t0 = 0, t = 0;
    double y0 = 1, y = 0;
    double h = 0.1;
    
    printf("Equation dx/dt = -x solution, when x(0) = 1 \n \nt    x       exp(-t) \n \n");
    
    t = t0;
    y = y0;
    
    for(int i = 0; i <= 10; i++)
    {
        printf("%3.1lf  %5.4lf  %5.4lf \n", t, y, exp(-t));
        step(&t, &y, h);
    }
}

void step(double *t, double *y, const double h)
{
    /* 4-th ORDER
     double k1 = f(*y);
    double k2 = f(*y + h/2*k1);
    double k3 = f(*y + h/2*k2);
    double k4 = f(*y + h*k3);
    *y = *y + h/6 * ( k1 + 2*k2 + 2*k3 + k4);
     */
    
    // 2-nd ORDER
    *y = *y + h*f( *y + 1/2.0*h*f(*y) );
    
    // INCREMENT t
    *t = *t + h;
}

double f(const double y)
{
    return -y;
}

Equation dx/dt = -x solution, when x(0) = 1 
 
t    x       exp(-t) 
 
0.0  1.0000  1.0000 
0.1  0.9050  0.9048 
0.2  0.8190  0.8187 
0.3  0.7412  0.7408 
0.4  0.6708  0.6703 
0.5  0.6071  0.6065 
0.6  0.5494  0.5488 
0.7  0.4972  0.4966 
0.8  0.4500  0.4493 
0.9  0.4072  0.4066 
1.0  0.3685  0.3679 
Program ended with exit code: 0