C++学习日记#4——改进欧拉(Euler)公式

在x，y的平面上，微分方程y'=f(x,y)的解y=y(x)称作它的积分曲线。积分曲线上的一点(x,y)的切线斜率等于函数f(x,y)值。如果按函数f(x,y)在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致。

在求解的方法中有：欧拉公式、梯形方法。梯形方法虽然提高了精度，但是其算法复杂，在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数f(x,y)的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测。——《数值分析 （第五版）李庆杨编著》

改进欧拉公式

先用欧拉公式求得一个初步的近似值，称之为预测值，预测值的精度可能很差，再用梯形公式将他矫正一次，校正结果称为校正值。建立的预测——校正系统称为改进欧拉公式。

                                     预测：

                                     校正：

为了便于编写程序通常表示为平均化的形式

在该日志中的实现程序是以改进欧拉公式平均形式为主函数进行求解

程序主要思路

三个主要模块：赋值模块、求解函数模块以及结果输出模块；

主要讲解求解函数模块

主函数for循环

与以往的for循环语句(前4篇学习日志)不同，这里的for循环以每次跨步后的x的值与最右端点比较 当x值大于最右端点的值时，结束循环。这样的好处：

1.每进行一次循环可以得到该循环的x值；

2.当可行域的区间长度不能整除步长时，可以计算出在可行域区间的x值以及y值；

调用函数

在编写程序时，发现yp与yc的求解公式相近，思考如何使用一个调用函数进行两次求解得出yp和yc。如果在求解yc时继续使用调用求解yp的函数，会出现原本是yn的值进行运算将会被yp替代，这是因为求解yc的公式只有三个参数量（yn、xn、h），而求解yp的公式有四个参数量（yn、xn+1、yp、h）。所以本程序中为了避免这种情况的发生，将调用函数统一为四个参数量的求解函数。当求解yp时，多出来的那个参数将赋值零（yn、xn、0、h），求解yc时，函数初始参数为（yp、xn+1、c、h）。。其中的C是调用函数前计算得出，C为yn减去yp的结果。

好处是：只有一个调用函数，在修改函数公式时，只要将h*f(x,y)的f(x,y)改为新的求解函数公式即可

double func_y(double x, double y, double c, double h)//定义平均化形式的改进欧拉公式yp 
{
    f_y = y + c + h * (pow(x,3)*pow(y,2) + 1);
//第一次调用c=0，
//结果yn + h * (pow(xn,3)*pow(yn,2) + 1);
//第二次调用c=yn-yp
//结果yp+(yn-yp) + h * (pow(xn+1,3)*pow(yp,2) + 1);
//即：yn + h * (pow(xn+1,3)*pow(yp,2) + 1)
    return f_y;
}

附上实现程序（Visual Studio 2019）

#include <iostream> 
#include <iomanip> //头文件——该程序中主要为了保留小数点后几位
#include <cmath>//头文件——可以使用常用的数学函数

using namespace std;
double x[100];//定义每次跨越一步的x值； 
double y[100];//定义每个x值对应的初值y 
double yp[100];//定义一个预测值yp
double C;//定义y与yp的差值，便于后面调用 
double func_y(double x, double y, double c, double h);
//定义平均化形式的改进欧拉公式(加入了一个c变量
//这个变量好处在于可以将一个调用函数只要改变参数就可以求出yp和yc) 
double f_y;//定义调用函数的返回值yp 
double a;//定义x可行域的下限 
double b;//定义x可行域的上限 
double h;//步长； 
int I;
//定义步数
//这主要是由于区间长度能否被步长整除是不可知的
//由于要保证每次跨步要在x的定义域内，步数是不可知的
//这里通过每次跨步后的x的值与最右端点比较 当x值大于最右端点的值，结束求解。 

int main()
{
    cout << "输入x的上限" << endl;
    cin >> b;
    cout << "输入x的下限" << endl;
    cin >> a;
    cout << "输入x0" << endl;
    cin >> x[0];
    cout << "输入y0" << endl;
    cin >> y[0];
    cout << "输入步长" << endl;
    cin >> h;//输入模块；之后为计算在可行域内，以h的步长可以最多行进多少次； 
    for (int i = 1; x[0] + i * h <= b; i++)
        //这里通过每次跨步后的x的值与最右端点比较 当x值大于最右端点的值，结束循环
    {
        x[i] = x[0] + i * h;
        func_y(x[i - 1], y[i - 1], 0, h);//第一次调用求出预测值yp 
        yp[i] = f_y;
        C = y[i - 1] - f_y;
        //这个变量的主要作用是为了统一调用函数的参数个数
        //可以将一个调用函数只要改变参数就可以求出yp和yc
        func_y(x[i], yp[i], C, h);//第二次调用求出yc 
        y[i] = (yp[i] + f_y) / 2;//yp与yc的平均值为矫正值 
        I = i;
    }
    cout << "**********输出结果**********" << endl;
    cout << setiosflags(ios::fixed) << setprecision(5);//输出结果保留5位小数 
    cout << "    Xn " << "    预测值  " << "  校正值" << endl;
    for (int i = 1; i <= I; i++)
    {
        cout <<" " << x[i] << " , " << yp[i] << " , " << y[i]  << endl;
    }
    return 0;
}

double func_y(double x, double y, double c, double h)//定义平均化形式的改进欧拉公式yp 
{
    f_y = y + c + h * (pow(x,3)*pow(y,2) + 1);
//修改新的求解函数时，只要改 (pow(x,3)*pow(y,2) + 1)括号里面的内容即可
    return f_y;
}

例题

输出结果

当可行域的区间长度不能整除步长时；

后话 

对前面的几篇日记进行了大致的回顾，发现都有一个共同点，定义的变量好像有点多，例如一些循环内部的临时变量以及很多数组。现在《C++ primer plus》自学到第四章，里面就讲解到了delete可以释放内存，结束第四章学习后，之后的程序可能有目的性的加入delete。还有就是有没有发现几篇日记定义变量时就只有两种类型double、int，这是一个现在对定义类型的一个盲区，现在不处理好着知识点的话，后面肯定会要花时间补上这一块知识；

在编写程序之前，其实也考虑过两个需求，1.可以自由输入简单公式；2.步长与区间长度的关系；

在实现第二个需求时，有考虑过if语句(主要判断能不能整除，能整除输出整数，不能整除输出商的整数部分，得出步数，之后就像往常的for循环求解)，在编写时，发现这样太麻烦了，想要判断效果，其实for语句中也有类似的效果，而且使用后发现不仅能求出步数，还可以将每个x赋值，一举多得；

关于第一个需求，其实没有实现。网络冲浪了一下，看到一些回答说不能实现的，好像只能通过调用组件来实现。之后我就没有继续研究下去了

现在《数值分析》已经结课了，之后会补上几篇书上的一些迭代公式，这些迭代公式大都是一些很简单的赋值、求解、输出三大块，其实没有应用到很多方面的知识。弄完之后就开始尝试做一下小项目。