**Question**

Consider an initial value ODE of the following form

As with all programs we start by thinking about what are the parameters and local variables in the problem. It is clear from the specification here that the parameters are a, b, n, as well as the function to be integrated itself although as we are not interested in writing a generic algorithm we can ignore the last one. Start with an empty program with libraries for input/output to screen/file and mathematical functions.

#include<iostream>#include<fstream>#include<cmath>usingnamespacestd; intmain() {return0; }

Compile and run the program to check you have everything set up ok.

Next shall declare the external parameters first, followed by the local parameters and then initialise any values at the start of the algorithm.

// parametersint n=10; double a=0.,b=1.,alpha=0.;// local variablesdouble h,x,w;// intialise valuesh=(b-a)/(double)(n); x=a; w=alpha;

Finally we can add in the rest of the algorithm in a couple of lines and output results to the screen at the same time.

// implement Euler's methodfor(int i=0;i<n;i++) { x = a + i*h;// update value of x to x_icout << x << " " << w << endl;// output x_i and w_i to screenw = w + h*(x*exp(3*x)-2.*w);// update w to w_{i+1}} cout << b << " " << w << endl;// output x_n and w_n to screen

The output from this program should give y(x=1)=2.7609. Now that the program is running we should check the value of the result for different values of n, and see if the results appear to converge to a constant value as the number of steps n increases.

Now we can move this code into its own function, so that we can be more flexible with obtaining results. If we use the external parameters as arguments to the function and then copy/paste the function in from ”main” it will look like:

doubleeulersMethod(int n,double a,double b,double alpha) {// local variablesdouble h,x,w;// intialise valuesh=(b-a)/(double)(n); x=a; w=alpha;// implement Euler's methodfor(int i=0;i<n;i++) { x = a + i*h;// update value of x to x_icout << x << " " << w << endl;// output x_i and w_i to screenw = w + h*(x*exp(3*x)-2.*w);// update w to w_{i+1}} cout << b << " " << w << endl;// output x_n and w_n to screenreturnw; }

and can be used inside ”main” like

cout << " y(b) ~ " <<eulersMethod(1000,0.,1.,0.) << endl;

If you want to output the results of x against w to a file for plotting then we can use the flexibility of the stream variable by including it as an argument to the function. We add this argument as ostream& output (must use a non constant reference) and then where cout appears in the function it can be replaced by output. The final code could look something like:

#include<iostream>#include<fstream>#include<cmath>usingnamespacestd; doubleeulersMethod(int n,double a,double b,double alpha,ostream& output) {// local variablesdouble h,x,w;// intialise valuesh=(b-a)/(double)(n); x=a; w=alpha;// implement Euler's methodfor(int i=0;i<n;i++) { x = a + i*h;// update value of x to x_ioutput << x << " , " << w << endl;// output x_i and w_i to screenw = w + h*(x*exp(3*x)-2.*w);// update w to w_{i+1}} output << b << " , " << w << endl;// output x_n and w_n to screenreturnw; } intmain() {// output to screeneulersMethod(1000,0.,1.,0.,cout);// output to fileofstreammyFileStream("test.csv");eulersMethod(1000,0.,1.,0.,myFileStream); myFileStream.close();return0; }

#### Euler Method: Convergence and Accuracy

Next we want to analyse the accuracy of the method. In this case we do have an analytic solution to compare against but often we don’t, so what can we do in that case? First, assume that if a method that converges at the rate c the following holds true

^{2}n steps where k is an integer. The ratio of differences between the estimates R, written as

To generate an empirical estimate for the convergence rate and confirm accuracy of the
method we need a slimmed down version of the previous function. Create a new project with the
eulerMethod function as follows

doubleeulersMethod(int n,double a,double b,double alpha) {// local variablesdouble h,x,w;// intialise valuesh=(b-a)/(double)(n); x=a; w=alpha;// implement Euler's methodfor(int i=0;i<n;i++) { x = a + i*h;// update value of x to x_iw = w + h*(x*exp(3.*x)-2.*w);// update w to w_{i+1}}returnw; }

Now we want to run a loop over different values of n to generate a table of results showing the convergence rate. To calculate the ratio R you will need to keep track of the difference between subsequent results. To do this you need to store the old version of value and difference outside the loop. The common way to denote these is valueOld and diffOld. Inside the loop store the value and difference using current value of n. The algorithm should look a bit like this

// evaulate the convergence of the method// using an empirical method, relies on there// being smooth convergence// ratio between number of steps at each stageint k=4;// store previous values and differencesdouble valueOld=1.,diffOld=1.;for(int i=1;i<=10;i++) { int n=pow(k,i);// calculate value with ndouble value =eulersMethod(n,0.,1.,0.);// and difference from last timedouble diff=value-valueOld;// output stage, steps, value, ratio R, convergence rate c and errorcout << i << " " << n << " " << value << " " << diffOld/diff <<endl;// store old valuesvalueOld = value; diffOld = diff; }

You might want to add in a bit of extra information to you outputs, the following code shows the convergence rate tending to 1 for the Euler method. Take this code and implement the midpoint rule or Runge-Kutta methods and check on the rate of convergence.

#include<iostream>#include<fstream>#include<cmath>usingnamespacestd; doubleeulersMethod(int n,double a,double b,double alpha) {// local variablesdouble h,x,w;// intialise valuesh=(b-a)/(double)(n); x=a; w=alpha;// implement Euler's methodfor(int i=0;i<n;i++) { x = a + i*h;// update value of x to x_iw = w + h*(x*exp(3.*x)-2.*w);// update w to w_{i+1}}returnw; } intmain() {// evaulate the convergence of the method// using an empirical method, relies on there// being smooth convergence// ratio between number of steps at each stageint k=4;// store previous values and differencesdouble valueOld=1.,diffOld=1.;for(int i=1;i<=10;i++) { int n=pow(k,i);// calculate value with ndouble value =eulersMethod(n,0.,1.,0.);// and difference from last timedouble diff=value-valueOld;// calculate Rdouble R=diffOld/diff;// can show that R = k^c where c is convergence ratedouble c=log(R)/log(k);// use c to estimate Adouble A=pow(n,c)/(1-pow(k,c))*diff;// output stage, steps, value, constant A, ratio R, convergence rate c and errorcout << i << " " << n << " " << value << " " << A << " " << R << " " << c << " " <<fabs(A/pow(n,c)) <<endl;// store old valuesvalueOld = value; diffOld = diff; }// output from this code// 1 4 2.09213 11.8537 0.915638 -0.063575 12.9459// 2 16 2.93241 -4.73586 1.29973 0.189106 2.80345// 3 64 3.14744 -4.41281 3.90769 0.983158 0.0739528// 4 256 3.20119 -4.58885 4.00079 1.00014 0.0179111// 5 1024 3.21462 -4.59139 4.00162 1.00029 0.00447473// 6 4096 3.21798 -4.58666 4.00049 1.00009 0.00111896// 7 16384 3.21882 -4.58472 4.00013 1.00002 0.000279766// 8 65536 3.21903 -4.58409 4.00003 1.00001 6.99432e-05// 9 262144 3.21908 -4.58391 4.00001 1 1.74859e-05// 10 1048576 3.21909 -4.58385 4 1 4.37148e-06return0; }

## No comments:

## Post a Comment