查看原文
其他

[C++] 中值法求函数值

2016-12-15 Y叔 biobabble

这个方法也可以叫二叉搜索,因为和使用二叉树进行搜索的过程是一样的,看值大小来决定是往右边走还是左边走,最后到了最接近的值,也就是答案了,返回。


思路很简单,随便搞两个起始数字(left, right),算出来函数值一个小于0,一个大于0, 那么再拿这两数的均值(M),再算一下函数值,如果>0,那么缩小范围到[left, M],反之则缩小到[M, right],如此反复,直到函数值小于某个tolerance,返回答案M,就是这么简单而已。


假设有以下的函数:

A * x + B * sqrt(x ^ 3) - C * exp(-x / 50) - D = 0

精度要求:0.0000001 = 1e-7


那么给定下面两组数字, ABCD:

input data: 2 0.59912051 0.64030348 263.33721367 387.92069617 15.68387514 1.26222280 695.23706506 698.72384731

我们将得到答案:

answer: 73.595368554162 41.899174957955


#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <cmath>

double f(double x, double A, double B, double C, double D);
double bs(double left, double right, double eps, double A, double B, double C, double D);

int main() {
  std::ifstream infile("data/id34.txt");
  std::string line;
  getline(infile, line);
  std::stringstream ss(line);
  int N;
  ss >> N;
  std::cout << std::setprecision(16); // show 16 digits
  double A, B, C, D;
  double eps = 0.0000001;
 
  for (int i=0; i<N; i++) {
    getline(infile, line);
    std::stringstream params(line);
    params >> A;
    params >> B;
    params >> C;
    params >> D;
    
    std::cout << bs(0.0, 100.0, eps, A, B, C, D) << " ";
    
  }
  std::cout << std::endl;
}

主程序就是在读数字,然后传给bs函数(binary search)求函数值,并打印结果。


首先定义函数,这个函数我们知道在0和100之间有解,这就给我们中值搜索带来极大便利。

double f(double x, double A, double B, double C, double D) {
  // A * x + B * sqrt(x ^ 3) - C * exp(-x / 50) - D = 0
  double res = A * x + B * pow(x, 1.5) - C * exp(-0.02*x) -D;
  return res;
}


然后这个中值求解,或者叫二叉搜索就很容易写:

double bs(double left, double right, double eps, double A, double B, double C, double D) {

  double M = (left+right)/2;
  double v = f(M, A, B, C, D);

  if (std::abs(v) <= eps) {
    return M;
  }
  if (v > 0) {
    right = M;
  } else {
    left = M;
  }
 
  bs(left, right, eps, A, B, C, D);
}

求中值,再求函数值,在精度范围内,则返回方程的解,否则继续寻找,无非是个递归。


这里有个小小的问题,cmath里的abs函数要用std::abs来调用,如果直接写abs会调用cstdlib的abs,只会返回int,小小的坑,花了我一点时间找出来。




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存