decay process with euler algorithm

[glao]

Hello , I am trying to solve some problems from the "Computer Simulation Methods" book.

It says

Write a class that solves the nuclear decay problem.Input the decay constant l from the control window.For l=1 and dt=0.01 ,compute the difference between the analytical result and the result of the euler algorithm for Nt/N0 at t=1 and t=2.Assume that time is measured in seconds.

(The above is for java ,but I want to do it in c++).

My question is how I will implement Euler algorithm.
I know ,that in motion for example , the Euler algorithm goes as :

void Euler(){
             x=x0+v0*dt;  
             v=v0-k*x0*dt;
             t=t+dt;
             error = fabs((x0-analyticPosition(1,0.01))/x0);
             v0=v;
             x0=x;

}

where x is the position ,v is velocity and xo,vo are initial values.

I can't understand in the decay situation where Nt/N0=exp(-l) how should I implement Euler.

Thanks!

[Moschops]

The simple, forward-Euler method works as follows:

1) Get current value of the variable you care about.
2) Get rate-of-change of the variable.
3) Next value is current_value + (rate_of_change * time_step)

Repeat until you've gone forwards enough in time.

So in this case, it looks like the time_step is 0.01, and the rate-of-change is (negative) one multiplied by the current_value.

Start with amount of stuff X at time t=0. So how much stuff is there at time t=0+0.01? X+ (-1*X*0.01)

[glao]

1) I care about Nt/N0 , so I will have an input for that variable.

2) I think that the problem wants Nt/N0 and not dN/dt as you say (if I am not wrong about you) ,so in this case how I wil compute / input the rate of change ?

3) So , I will start with an initial amount of Nt/N0 ?(if the above is correct).

Also , because the problem wants Nt/N0 ,it is not a negative value.Negative is dN/dt , right?

Thank you.

[Moschops]

1) I care about Nt/N0 , so I will have an input for that variable.

Nt/N0 is what you're calculating. It's the amount of stuff at time t, divided by the amount of stuff at time zero (i.e. at the start). You care about the amount of stuff.

how I wil compute / input the rate of change ?

As I said, the rate of change is negative one multiplied by the current amount of stuff. The question says I=1. What's the rate of change at the start? -1 * N0. What's the rate of chnage 0.01 seconds later? -1 * N(t=0.01) ; minus one times the amount of stuff at time 0.01

3) So , I will start with an initial amount of Nt/N0 ?(if the above is correct).

No. That's what you're calculating. Nt is the amount of stuff left at time t. N0 is the amount of stuff at the start. You will start with an inital amount of N, called N0.

[glao]

Ok, so I will have sth like:

dt=0.01;
N0=0;

void Euler(){
             N=N0-1*N0*dt;
             t=t+dt;
             error = fabs(N0-analyticalDecay(1,0.01))/N0);
             N0=N;


}

double analyticalDecay(int l,int dt){

return exp(-l*dt);


}

I am not sure about the

N=N0-1*N0*dt;

I am confused because the analytical is N=N0exp(-l*t)

Thanks

[mike_2000_17]

Do you know calculus?

If you have a first-order differential equation like this:

dN/dt = -l * N(t)

You get the following analytical solution:

N(t) = N(t0) * exp(-l*t)

Here's how you get it:
1) Separate the variables:

dN / N(t) = -l * dt

2) Take the integral of both sides:

ln( N(t) ) = -l * t + C

3) Isolate the dependent variable:

N(t) = exp( -l * t + C ) = K * exp( -l * t )

4) Apply the initial condition N(0) = N0:

N(0) = K * exp( 0 ) = K = N0

N(t) = N0 * exp( -l * t )

And that's it. This is the most elementary form of ODE (Ordinary Differential Equation), a particular kind ODE that you usually learn in high-school calculus. Methods like Euler integration are methods used to solve these ODE problems numerically, because most of them cannot be solved analytically, except for the simplest ones like the above example.

And yes, your solution:

N=N0-1*N0*dt;

is correct.

[glao]

Οκ, thank you both very much!

[glao]

I have a problem.The numerical and analytical solution is exactly the same no matter the number of iterations.

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <iomanip>

using namespace std;


class Decay {
private:
    double l,dt,Nt,N0,t;
    double  error=10.0;//initial value to error in order to compare with tolerance
    double  eps=1e-6; //my error tolerance
public:
    Decay(){};
    void read();
    void Euler(){
             Nt=N0-1*N0*dt;
             //Nt=N0*exp(-l*dt);
             t=t+dt;
             error = fabs((N0-analyticalDecay(N0,l,dt))/N0);
             N0=Nt;
}

double analyticalDecay(double N0,int l,int dt){

return N0*exp(-l*dt);


}


 void results(){
            cout <<"\nNumerical results: \n";
            cout <<"Nt = " <<Nt <<"\n";
            cout<<"\nAnalytical Decay: \n"<<"Nt = "<<analyticalDecay(N0,l,dt);
            cout <<"\nThe error in amount N : "<<error*100 <<"%"<<endl;

        }


};

int main()
{
Decay mydecay;

mydecay.read();

for (int i=0;i<10000;i++)
mydecay.Euler();


cout <<"\nResults : "<<endl;
mydecay.results();



    return 0;
}


void Decay::read(){

   cout <<"\nEnter constant l in seconds : "<<endl;
   cin >>l;
   cout <<"\nEnter dt step in seconds : "<<endl;
   cin >>dt;
   cout <<"\nEnter N0 : "<<endl;
   cin >>N0;

}

Where am I wrong?

Thanks!

[Moschops]

Shouldn't you be calcuating your analytical result with the original value of N0, rather than the value it is when your numerical method has finished with it?

[glao]

Ok, so I did:

class Decay {
private:
    double l,dt,Nt,N0,t;
    double initialN0=10; //I added this as the original amount
    ...

    double analyticalDecay(double initialN0,int l,int dt){

return initialN0*exp(-l*dt);

}

void results(){
            ...
            cout<<"\nAnalytical Decay: \n"<<"Nt = "<<analyticalDecay(initialN0,l,dt);

but it gives me always analytical result the value 10.(the initialN0)

[glao]

Any ideas here?Thanks!

[glao]

Hello , I will appreciate any help, thanks~

[Moschops]

Shouldn't the initialN0 value be whatever the user types in, instead of always 10?

Also, dt is your timestep. Shouldn't the analytical decay calculation be using the whole time, t?

Why do you calculate error ten thousand times?

[glao]

Hello and thank for the help.

I tried this (I use initialN0=700 and I input No 700 just for testing ) , but it returns me N=0 for the analytical.
Also , I am not sure about error as you said.I can't figure how to use it.

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <iomanip>

using namespace std;


class Decay {
private:
    double initialN0=700;
    double l,dt,Nt,N0,t;
    double  error=10.0;//initial value to error in order to compare with tolerance
    double  eps=1e-6; //my error tolerance
public:
    Decay(){};
    void read();
    void Euler(){
             Nt=N0-1*N0*dt;
             //Nt=N0*exp(-l*dt);
             t=t+dt;
             error = fabs((N0-analyticalDecay(initialN0,l,dt*100000))/N0);
             N0=Nt;
}


double analyticalDecay(double initialN0,int l,int t){
t=dt*100000;
return initialN0*exp(-l*t);


}


 void results(){
            cout <<"\nNumerical results: \n";
            cout <<"Nt = " <<Nt <<"\n";
            cout<<"\nAnalytical Decay: \n"<<"Nt = "<<analyticalDecay(initialN0,l,dt*100000);
            cout <<"\nThe error in amount N : "<<error*100 <<"%"<<endl;

        }


};

int main()
{
Decay mydecay;

mydecay.read();

for (int i=0;i<100000;i++)
mydecay.Euler();


cout <<"\nResults : "<<endl;
mydecay.results();



    return 0;
}


void Decay::read(){

   cout <<"\nEnter constant l in seconds : "<<endl;
   cin >>l;
   cout <<"\nEnter dt step in seconds : "<<endl;
   cin >>dt;
   cout <<"\nEnter N0 : "<<endl;
   cin >>N0;

}

[Moschops]

Well now you appear to be using the value dt*100000*100000 in your analytical calculation, because you multipy by 100000 before you pass it in, and then again afterwards. This seems wrong.

Set t to zero at the start.

Error is meant to be the difference between your analytical calcualtion, and what your Euler method came up with. So here is an example of how to calculate the error:

error = analyical_answer - euler_answer;

You need to get your final value for N from the Euler, and the correct answer from the analytical function, and subtract them. This calculation needs to be done once.

[glao]

Ok, I corrected:

double analyticalDecay(double initialN0,int l,int t){
//t=dt*100000;  I removed this
return initialN0*exp(-l*t);


}

I set double l,dt,Nt,N0,t=0; (t=0)

and

void results(){

            error=fabs(analyticalDecay(initialN0,l,t)-N0);
            cout <<"\nNumerical results: \n";
            cout <<"Nt = " <<Nt <<"\n";
            cout<<"\nAnalytical Decay: \n"<<"Nt = "<<analyticalDecay(initialN0,l,t);
            cout <<"\nThe error in amount N : "<<error*100 <<"%"<<endl;

        }

but it gives me always Nt=0 (analytical result)

I can't figure it!

[Moschops]

but it gives me always Nt=0 (analytical result)

That's the correct answer after 1000 seconds; it's pretty much all decayed. You're calculating to 1000 seconds. You should be using 100 steps to get to t=1, and another 100 to get to t=2. Why are you using 100000?

Why are you using int values of t and l in your analytical function, but double values of them elsewhere? If they're not whole numbers, you're going to be losing information when you convert them to an int.

I see that you set initialN0 to 700 and then never change it, but you still use it in your calculattions. This makes no sense. What also makes no sense is asking the user to give an N0 amount. You don't care what value N0 is so long as it's not zero at the start.

[glao]

Ok!I totally missed that about the time!Of course they all decay!

I also changed the values to double (I have forgotten that).

Now , I set initialN0=700 for the analytical result and N0 for the user (in which I enter again 700 to test) because we said before that the analytical result should contain the original value of N0.
The user enters N0 because I must know the value at the beginnig.
Note ,that I have initialN0 for analytical and N0for user.

Because ,if I use the double analyticalDecay(double N0,.. instead of double analyticalDecay(double initialN0,...{
it will compute N0 as the last value from euler steps as you said before,right?

[Moschops]

Yes. Both methods need to start with the same initial value of stuff. it just doesn't matter what that initial value is, as long as it isn't zero.

[glao]

So , I think this code is ok.

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <iomanip>

using namespace std;


class Decay {
private:
    double initialN0=700;
    double l,dt,Nt,N0,t=0;
    double  error=10.0;//initial value to error in order to compare with tolerance
    double  eps=1e-6; //my error tolerance
public:
    Decay(){};
    void read();
    void Euler(){
             Nt=N0-1*N0*dt;
             //Nt=N0*exp(-l*dt);
             t=t+dt;

             N0=Nt;
}


double analyticalDecay(double initialN0,double l,double t){
//t=dt*100000;
return initialN0*exp(-l*t);


}


 void results(){
            error = fabs((N0-analyticalDecay(initialN0,l,t))/N0);
            cout <<"\nNumerical results: \n";
            cout <<"Nt = " <<Nt <<"\n";
            cout<<"\nAnalytical Decay: \n"<<"Nt = "<<analyticalDecay(initialN0,l,t);
            cout <<"\nThe error in amount N : "<<error*100 <<"%"<<endl;

        }


};

int main()
{
Decay mydecay;

mydecay.read();

for (int i=0;i<100;i++)
mydecay.Euler();


cout <<"\nResults : "<<endl;
mydecay.results();



    return 0;
}


void Decay::read(){

   cout <<"\nEnter constant l in seconds : "<<endl;
   cin >>l;
   cout <<"\nEnter dt step in seconds : "<<endl;
   cin >>dt;
   cout <<"\nEnter N0 : "<<endl;
   cin >>N0;

}

Thanks for your help!

[Moschops]

So if the user enters a value of 10 for N0, what value of initialN0 do you use in your analytic function?

[glao]

I know I must have the same value for these.
I did this :

class Decay {
private:
    double initialN0=0;
    ......

void Decay::read(){

   cout <<"\nEnter constant l in seconds : "<<endl;
   cin >>l;
   cout <<"\nEnter dt step in seconds : "<<endl;
   cin >>dt;
   cout <<"\nEnter N0 : "<<endl;
   cin >>N0;
   initialN0=N0;  // added this line

}

Is there a better way to do this?

Thank you.

[glao]

Hello ,if you have any ideas I want to hear.

Thanks!