关于数学分析程序设计作业快速入门
数分骚骚函数作业快速入门
前言
我也不知道自己写的对不对， 也没有数据来测， 所以我欢迎一切质疑和讨论。
本人水平也有限， 很多地方讲的不到位， 所以说可以直接来找到我来讨论！
本篇由C语言编写 ， 因为我上的是C语言程序设计 （
（毕竟有一个东西叫知识的诅咒qwq）
----by ska_0x08 2026 4 24
前置
什么是 typedef double (*func_t)(double x);
第一个 double 代表的是该函数返回的是一个double 数值
第二个 double 代表的是传入该函数的值是一个double 类型
我们称 “传入一个double 传出一个double” 这种类型的函数为 *func_t
在这个语句之后， func_t 就是一种声明
例如
func_t f;
意思就是 我定义了一个 func_t 类型的函数, 他的名字叫 \(f\) , 我在使用它的时候要传入一个double , 它吐出一个 double 出来。
下面是一个具体的例子
#include <stdio.h>
typedef double (*func_t) (double x);
double cube (double x){
return x * x * x;
}
int main (){
func_t p = cube;
printf("%f\n", p(3));
}
事实上 类似于
double cube (double x){
return x * x * x;
}
的部分应该在老师那里， 也就是说我们并不知道具体函数是什么， 但是我们可以使用 \(f(x_0)\) 这样的调用来获取在 \(x_0\) 处的值。 事实上， 对于本次的编程任务， 我们也不需要具体的函数， 只需要知道它在几点处的取值即可。
正片开始
Part.1 \(\space\) composite_trapezoid
具体要求
* 复化梯形公式
* @param f 被积函数
* @param a 积分下限
* @param b 积分上限
* @param m 等分数（子区间个数）
* @return 积分近似值
分析
\(h = \frac{b - a}{m}\)
\(T_m^{(1) }= \frac{h}{2}\sum\limits_{i = 1}^{m}[f(x_{i - 1}) + f(x_i)]\)
\(x_i = a + ih\)
我们要求解 \(Ans\) 就需要 \(x_i\) , \(x_i\) 需要 \(a, h\) 而 \(h\) 只需要 \(b, a, m\) 即可求解。
为了增加代码可读性 , 我们用 数组来存储, \(x[i]\) 就是 \(x_i\) 的值
My coding
double composite_trapezoid(func_t f, double a, double b, int m) {
double x[114514191];
double h = (b - a) / m;
for (int i = 0; i <= m; i ++) //这里一定从0开始
x[i] = a + i * h;
double ans = 0;
for (int i = 1; i <= m; i ++)
ans += (f(x[i - 1]) + f(x[i]));
ans = ans * h / 2.0;
return ans;
}
Part.2 \(\space\) composite_simpson
具体要求
* 复化 Simpson 公式
* @param f 被积函数
* @param a 积分下限
* @param b 积分上限
* @param m 等分数（每段再 2 等分做 Simpson）
* @return 积分近似值
分析
由等式
\(T_m ^{(2)} = \frac{4T_{2m}^{(1)} - T_M^{(1)}}{4 - 1}\)
也就是说我们直接调用
第一题的答案即可。
double composite_simpson(func_t f, double a, double b, int m) {
return (4.0 * composite_trapezoid(f, a, b, 2 * m) - composite_trapezoid(f, a, b, m)) / 3.0;
}
我不知道这样子能否通过本题， 所以我们还是从头好好写一遍吧！
\(T_m^{(2)} = \frac{h}{6} [f(a) + f(b) + 2 \sum\limits_{i = 1}^{m - 1}f(x_i) + 4 \sum\limits _{i = 1}^{m} f(x_{i - \frac{1}{2}})]\)
其中根据定义我们知道
\(x_{i - \frac{1}{2}} = \frac{x_{i - 1} + x_i}{2}\)
那么也就可以根据公式写出来啦
My coding
double composite_simpson(func_t f, double a, double b, int m) {
double x[114514191];
double h = (b - a) / m;
for (int i = 0; i <= m; i ++) //这里一定从0开始
x[i] = a + i * h;
double ans = 0;
ans += f(a);
ans += f(b);
for (int i = 1; i <= m - 1; i ++)
ans += 2.0 * f(x[i]);
for (int i = 1; i <= m; i ++){
double x_ = (x[i - 1] + x[i]) / 2.0;
//这里x_ 对应于 下表为 i-(1/2) 的x的值
ans += 4.0 * f(x_);
}
ans = ans * h / 6.0;
return ans;
}
小注 : 前两种方法大抵可以不用数组把 \(x_i\) 存储下来， 我在此处这么写只是为了提高可读性。
Part.3 \(\space\) romberg
具体要求
* Romberg 方法
* @param f 被积函数
* @param a 积分下限
* @param b 积分上限
* @param tol 收敛精度（对角线相邻元素之差的绝对值）
* @return 积分近似值
分析
公式
\(T_m ^ {(k + 1)} = \frac{4^kT_{2m}^{(k)} - T_{m} ^ {k}}{4^k - 1}\)
这个实现起来有一定难度， 我们来逐步分解。
根据任务需求， 我们只有在 \(\vert T_m^{(k)} - T_m^{(k-1)}\vert < tol\) 的时候才计算结束。
我们考虑这么构建代码
- 计算 \(T_m^{(k)}\)
- 若 \(\vert T_{m} ^ {(k)} - T_m ^ {(k - 1)}\vert < tol\) 结束循环， 返回 \(T_m ^{(k)}\)
- 否则 \(k = k + 1\) 进入下一次循环。
下面的重头戏在于如何计算 \(T_m ^{(k)}\)
为了避免重复计算 我们将之前计算的值存储起来
我们记录二维数组 \(T\) 。 根据书上 \(P_{300}\) 公式我们会发现我们只会计算形如 \(T_{2^p m} ^{(k)}\) 这样的式子， 所以我们记录 \(T[k][p] = T_{2^pm}^{(k)}\)
根据书上的指示， 在求解 \(T_m^{(k)}\) 的时候， 我们需要 \(T_{m} ^{(k - 1)}\) 以及 \(T_{2m}^{(k - 1)}\)
注意 ！ 这个时候 \(T_m ^ {(k - 1)}\) 我们之前已经求解过了！ 也就是说问题转换成了 求解 \(T_{2m} ^{(k - 1)}\)
而 \(T_{2m} ^ {(k - 1)}\) 需要 \(T_{2m}^{(k - 2)}\) 和 \(T_{4m} ^ {(k - 2)}\) , 注意 这个时候 \(T_{2m} ^ {(k - 2)}\) 已经是求解过的了。
反复上述过程我们发现 最终我们只需要求解 \(T_{2^{(k - 1)}m}^ {(1)}\) 即可
根据 \(T_{2^{k - 1}m}^ {(1)}\)和 \(T_{2^{k - 2}m}^ {(1)}\) 求解 \(T_{2^{k - 2}m}^ {(2)}\)
根据 \(T_{2^{k - 2}m}^ {(2)}\)和 \(T_{2^{k - 3}m}^ {(2)}\) 求解 \(T_{2^{k - 3}m}^ {(3)}\)
......
根据 \(T_{2^{k - i}m}^ {(i)}\)和 \(T_{2^{k - i - 1}m}^ {(i)}\) 求解 \(T_{2^{k - (i + 1)}m}^ {(i + 1)}\) （这个其实是我们编写程序的关键）
......
根据 \(T_{2^{1}m}^ {(k - 1)}\)和 \(T_{2^{0}m}^ {(k - 1)}\) 求解 \(T_{2^{0}m}^ {(k)}\)
直接循环处理即可
My Coding
double romberg(func_t f, double a, double b, double tol) {
double T[1141][1141];
T[0][0] = 0;
int k = 1;
int m = 91;
double pow2 = 1.0;
while (1){
T[1][k - 1] = composite_trapezoid(f, a, b, pow2 * m);
double pow4 = 4.0;
for (int i = 1; i <= k - 1; i ++){
T[i + 1][k - i - 1] = ((pow4 * T[i][k - i] - T[i][k - i - 1]) / (pow4 - 1));
pow4 *= 4.0;
}
if (fabs(T[k][0] - T[k - 1][0]) < tol) break;
k ++;
pow2 *= 2.0;
}
return T[k][0];
}
当然这样做对我来说还是不完美， 仔细观察过程， 我们会发现有很多空间被浪费掉了， 于是我们有一个节省空间的写法。
缺点就是不直观， 看得明白就看看吧（感兴趣的可以直接来问我ovo）也不是很难的 trick
double romberg(func_t f, double a, double b, double tol) {
double T[11451419];
T[0] = 0;
int k = 1;
int m = 91;
double pow2 = 1.0;
while (1){
double pre = composite_trapezoid(f, a, b, pow2 * m);
double tmp;
double pow4 = 4.0;
for (int i = 1; i <= k - 1; i ++){
tmp = pre;
pre = ((pow4 * pre - T[i]) / (pow4 - 1));
T[i] = tmp;
}
T[k] = pre;
if (fabs(T[k] - T[k - 1]) < tol) break;
k ++;
pow2 *= 2.0;
}
return T[k];
}
Part.4 \(\space\) gauss_legendre
具体要求
* Gauss-Legendre 求积公式
* @param f 被积函数
* @param a 积分下限
* @param b 积分上限
* @param n Gauss 点个数
* @return 积分近似值
分析
本质是
\(ans = \frac{b - a}{2}\sum\limits_{i = 1} ^ {n + 1} a^*_i f(\frac{b-a}{2}x^*_i + \frac{a+b}{2})\)
其中当 \(n\) 固定之后， 其对应的 \(x^*_i\) 和 \(a^*_i\) 都是固定的。
具体见书 \(P_{302}\) 页。
在实现的时候， 我们通过打表把 \(n\) 所确定的 \(x^*_i\) 和 \(a^*_i\) 数列预处理出来， 然后就直接根据上式求值即可。
My Coding
double gauss_legendre(func_t f, double a, double b, int n) {
double x_1[3] = {0, 0.5773502692, -0.5773502692}, a_1[3] = {0, 1.0, 1.0};
double x_2[4] = {0, 0.7745966692, -0.7745966692, 0}, a_2[4] = {0, 0.5555555556, 0.5555555556, 0.8888888889};
double x_3[5] = {0, 0.3399810436, -0.3399810436, 0.8611363116, -0.8611363116},
a_3[5] = {0, 0.6521451549, 0.6521451549, 0.3478548451, 0.3478548451};
double x_4[6] = {0, 0.5384693101, -0.5384693101, 0.9061798459, -0.9061798459, 0},
a_4[6] = {0, 0.4786286705, 0.4786286705, 0.2369268851, 0.2369268851, 0.5688888889};
double ans = 0;
double tmp1, tmp2;
tmp1 = (b - a) / 2.0;
tmp2 = (b + a) / 2.0;
if (n == 1){
for (int i = 1; i <= n + 1; i ++){
x_ = x_1[i] * tmp1 + tmp2;
ans += f(x_) * a_1[i];
}
}
if (n == 2){
for (int i = 1; i <= n + 1; i ++){
x_ = x_2[i] * tmp1 + tmp2;
ans += f(x_) * a_2[i];
}
}
if (n == 3){
for (int i = 1; i <= n + 1; i ++){
x_ = x_3[i] * tmp1 + tmp2;
ans += f(x_) * a_3[i];
}
}
if (n == 4){
for (int i = 1; i <= n + 1; i ++){
x_ = x_4[i] * tmp1 + tmp2;
ans += f(x_) * a_4[i];
}
}
ans *= tmp1;
return ans;
}
尾
相信大家已经跃跃欲试了， 快点来试试吧！
浙公网安备 33010602011771号