尤拉方法 常微分

C++ Coding - Euler's Method






Question

Consider an initial value ODE of the following form





write a function to return the value of y at x=b given





For this example we are going to implement the Euler method as given by







where w at the nth step gives an estimate for the value of y at x=b.

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>
using namespace std;

double eulersMethod(int n,double a,double b,double alpha)
{
  // local variables
  double h,x,w;
  // intialise values
  h=(b-a)/(double)(n);
  x=a;
  w=alpha;
  // implement Euler's method
  for(int i=0;i<n;i++)
  {
    x = a + i*h; // update value of x to x_i
    w = w + h*(x*exp(3.*x)-2.*w); // update w to w_{i+1}
  }
  return w;
}

int main()
{
  
  // evaulate the convergence of the method
  // using an empirical method, relies on there
  // being smooth convergence
  // ratio between number of steps at each stage
  int k=4;
  
  // store previous values and differences
  double valueOld=1.,diffOld=1.;
  for(int i=1;i<=10;i++)
  {
    int n=pow(k,i);
    // calculate value with n
    double value = eulersMethod(n,0.,1.,0.);
    // and difference from last time
    double diff=value-valueOld;
    // calculate R
    double R=diffOld/diff;
    // can show that R = k^c where c is convergence rate
    double c=log(R)/log(k);
    // use c to estimate A
    double A=pow(n,c)/(1-pow(k,c))*diff;
    // output stage, steps, value, constant A, ratio R, convergence rate c and error
    cout << i << " " << n << " " << value << " " << A << " " << R << " " << c << " " << fabs(A/pow(n,c)) <<endl;
    
    // store old values
    valueOld = 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-06 
  return 0;
}

輸出結果


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-06