用C语言实现定积分计算(包括无穷积分/可自定义精度)

news2024/9/21 2:32:38

关于严谨性的声明:

在用C语言进行定积分的计算之前,我需要声明以下几点:

一、我们所进行定积分计算的函数都是应当是黎曼可积的,这保证了我们即使均匀地分割区间也保证了积分的收敛性。

二、我们同时还应该认识到,鉴于某些函数不一定是在每一点上都是可导的,因此我们将不使用所谓的利用到导数来提高计算精度的这种算法。

三、考虑到某些函数并不是连续的,如果涉及到这些计算,那么我们计算到的结果将会是一个瑕积分,如果这个瑕积分是收敛的,那么我们得到的就是收敛值。如果不收敛,这个瑕积分的值并不是积分主值,因为我们没有办法保证瑕点附近的两个点的趋近于瑕点的趋势是一样的。

四、如果你计算的是一个无穷积分,如果这个无穷积分是收敛的,那么我们得到的就是收敛值。如果不收敛,那么我们得到的值将会是它的积分主值。


另外,考虑到C语言精度的限制,我们无需写一个专门的程序用来计算瑕积分,因为当我们用C语言计算这个瑕积分的时候,就相当于我们采用了两边趋近瑕点取极限的办法来计算瑕积分,当然了,由于我们采用了均匀分割的方法,所以当我们两边趋近的时候并不能使得其以任意方式趋近瑕点。除此之外,你应当保证瑕点不在积分区间的边界上。对于无穷积分,不论我们怎么切割区间都不能使得划分的区间足够小,除非我们划分区间的时候将他划为了无穷份,当然了这个划分的份数应当是无穷积分区间大小的高阶无穷大。事实上,这是做不到的,因此我们将采用函数变换的方法将无穷积分转化为有限积分。

最后,鉴于C语言的精度限制,如果我们只是计算性质较为良好的函数,那么我们将会获得非常好的精度,甚至是仅仅只是可测的函数,也是可以计算的。但如果我们计算的是有着许多无穷间断点,或有着许多瑕点的函数,或者是震荡很厉害的函数,那么它的积分值误差将会较大。最后我还想声明一点,那就是我们不应该用该程序去计算本身就是发散的积分因为你甚至可能得到的仅仅只是个有限值。


1. 有限积分

根据我们在之前所说的,我们可以轻易得到:

4876ac26b780492fa8be38d3a18643af.png

 其中:

af0e38abe5074924a2496a1f9e8b34c5.png

 因此,我们不难写出以下代码:

//定义一个C语言函数,用于进行微积分的数值运算
long double integrate(long double a, long double b, long int n)
{
    long double s = 0;
    long double dx = (b - a) / n;
    for (long int i = 1; i < n; i++)
    {
        s += f(a + i * dx) * dx;
    }
    return s;
}

其中函数 long double integrate(long double a, long double b, long int n) 就是用来计算定积分的函数。值得注意的是,在调用该函数之前,我们应当提前定义被积函数f(x)。其中a,b,n分别为积分下限,积分上限和区间分割数。

以上是矩形法,数值计算的精度并不是非常高,我们改进一下,使用梯形法,得到:

2b8d53cf63924c3abcfe887ddc7511f3.png

  其中:

af0e38abe5074924a2496a1f9e8b34c5.png

  因此,我们不难写出以下代码:

long double integral(long double x0, long double xt, long long n)
{
    long double sum = 0;
    long double dx = (xt - x0) / n;

    for (long long k = 0; k < n; k++)
    {
        sum += ((f(x0 + k * dx) + f(x0 + (k + 1) * dx)) / 2) * dx;
    }

    return sum;
}

相较于矩形法,梯形法不论是精度还是速度都有着不小的提升。


2. 无穷积分

我们已经讨论过无穷积分采用直接计算的困难的地方,因此我们将采用函数变换的方法。

我们很容易注意到:

fc7ab638fa94417cbfd96d1a8366f04c.png

 

4b9f1397223e403983affb661e6d0c29.png

这两个式子是等价的,因此我们可以将所有的无穷积分转化为无限积分。只需要利用积分区间的可叠加性即可。因此不多叙述。


 

有了相关的理论铺垫后,我们添加一点点细节。以下直接给出关于定积分的计算代码:(详细信息在注释里)

#include <stdio.h>  
#include <math.h>  

#define inf INFINITY

// 定义用于被积分的函数  
long double f(long double x) 
{
    return exp(x);
}

// 定义被积函数的特殊形式用于计算无穷积分  
long double g(long double x) 
{
    long double y = 1.0 / x;
    long double result = f(y) * (1.0 / (x * x));
    return result;
}

// 一个用于计算定积分的C语言子程序
//值得注意的是我们还将数学函数作为参数传给了这个函数,这可以让我们有多个数学函数时,可以指定数学函数计算  
long double integral(long double x0, long double xt, long long n, long double (*f)(long double)) 
{
    //对于黎曼可积的函数,这样的积分值为零,直接返回。
    if ((x0 == xt)&&((x0 != INFINITY && x0 != -INFINITY)&&(xt != INFINITY && xt != -INFINITY)))
        return 0;
    else
    {
        long double sum = 0;
        long double dx = (xt - x0) / n; 

        for (long long k = 0; k < n; k++)
        {
            sum += ((f(x0 + k * dx) + f(x0 + (k + 1) * dx)) / 2) * dx;
        }

        return sum;
    }
}

// 定义绝对值函数  
static long double lfabs(long double x) 
{
    if (x >= 0)
        return x;
    else
        return -x;
}

//用于计算指定精度的定积分的函数
//值得注意的是我们还将数学函数作为参数传给了这个函数,这可以让我们有多个数学函数时,可以指定数学函数计算  
long double definite_integral(long double x0, long double xt, long double esp, long double (*f)(long double)) 
{
    //如果积分限不是无穷则进行经典数值积分计算
    if ((x0 != INFINITY && x0 != -INFINITY) && (xt != INFINITY && xt != -INFINITY)) 
    {
        long long i = 10;  
        long double d = 0;
        long double s = 0;

        do {
            long double s1 = integral(x0, xt, i, f);
            long double s2 = integral(x0, xt, 10 * i, f);

            i *= 10;
            s = s2;
            d = lfabs(s1 - s2);
        } while (d > esp);
        return s;
    }

    //如果是无穷积分,则采用特殊方法进行数值计算
    else 
    {
        //积分下限有限大小
        //值得注意的是,当我们使用g函数进行积分的时候,有一个积分限迎丹是0,然而事实上我们以一个很小的小数来替代
        if (x0 != INFINITY && x0 != -INFINITY)
        {
            if (xt == INFINITY)
            {
                long double s1 = definite_integral(x0, 1, esp, f);
                long double s2 = definite_integral(esp/1000000.0, 1, esp, g);
                return (s1 + s2);
            }
            else if (xt == -INFINITY)
            {
                long double s1 = definite_integral(-1, x0, esp, f);
                long double s2 = definite_integral(-1, -esp/1000000.0, esp, g);
                return -(s1 + s2);
            }
        }

        //积分上限有限大小
        else if (xt != INFINITY && xt != -INFINITY)
        {
            if (x0 == INFINITY)
            {
                long double s1 = definite_integral(xt, 1, esp, f);
                long double s2 = definite_integral(esp/1000000.0, 1, esp, g);
                return -(s1 + s2);
            }
            else if (x0 == -INFINITY)
            {
                long double s1 = definite_integral(-1, xt, esp, f);
                long double s2 = definite_integral(-1, -esp/1000000.0, esp, g);
                return (s1 + s2);
            }
        }

        //积分上下限都为无穷
        else if (x0 == -INFINITY && xt == INFINITY)
        {
            long double s1 = definite_integral(-1, 1, esp, f);
            long double s2 = definite_integral(-1, -esp/1000000.0, esp, g);
            long double s3 = definite_integral(esp/1000000.0, 1, esp, g);
            return (s1 + s2 + s3);
        }
        else if (x0 == INFINITY && xt == -INFINITY)
        {
            long double s1 = definite_integral(-1, 1, esp, f);
            long double s2 = definite_integral(-1, -esp/1000000.0, esp, g);
            long double s3 = definite_integral(esp/1000000.0, 1, esp, g);
            return -(s1 + s2 + s3);
        }

        //如果积分是从正无穷积分到正无穷或从负无穷积分到负无穷则需要另外判断
        else
        {
            if (x0 == INFINITY && xt == INFINITY)
            {
                if (f(INFINITY) == 0)
                    return 0;
                else
                    return NAN;
            }
            else if (x0 == -INFINITY && xt == -INFINITY)
            {
                if (f(-INFINITY) == 0)
                    return 0;
                else
                    return NAN;
            }
        }
        return NAN;
    }
}

int main() 
{
    long double x0, xt, esp;
    printf("如果要输入无穷大或无穷小,请输入inf或-inf\n");
    printf("请输入积分区域(x0,xt),和最大容忍误差esp:>>");
    scanf("%lf %lf %lf", &x0, &xt, &esp);

    long double result = definite_integral(x0, xt, esp, &f);
    printf("%.16lf\n", result);

    return 0;
}

在计算指定定积分精度的C函数中,对于无穷积分我们通过转化为有限积分后再调用自身的方式有效地减伤了大量大代码量。这也是为什么要把我们定义的数学函数作为参数传给函数的原因。

如下:

c94ba1cf189e45c69bc20fd96075e074.png

 


代码的运行:

3ffd0c7a889343b6a048668e63298a06.png

 我们指定数学函数为exp(x),可以看到我们很好地计算出来该无穷积分。

 

 

 

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/856100.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

一文详解 DolphinDB SQL 标准化

为了提升用户体验&#xff0c;降低用户学习成本和脚本迁移复杂度&#xff0c;自 1.30.17 / 2.00.5 版本开始&#xff0c;DolphinDB 逐步支持了标准化 SQL 的书写方法&#xff1b;并于 1.30.22 / 2.00.10 版本起&#xff0c;对标准 SQL 的常用语法和关键字实现了兼容。 1. 与标…

06-5_Qt 5.9 C++开发指南_Splash 与登录窗口(MouseEvent鼠标事件;注册表;加密)

文章目录 1. 实例功能概述2. 对话框界面设计和类定义3. QDIgLogin 类功能实现3.1 构造函数里的初始化3.2 应用程序设置的存储3.3 字符串加密3.4 用户名和密码输入判断3.5 窗口拖动功能的实现 4. Splash 登录窗口的使用5. 框架及源码5.1 可视化UI设计5.2 qdlglogin.h5.3 qdlglog…

使用vue模拟通讯录列表,对中文名拼音首字母提取并排序

一个功能需求,做一个类似联系人列表的功能,点击名称获取对应的id,样式简陋,只是一个模板,原来是uniapp项目,根据需要改成了vue,需要的自行设计css&#xff08;万是有一个mo的音&#xff09; 流程 获取数据提取首个字的拼音的首个字母排序并分组 上代码&#xff1a; <temp…

zadig安装驱动潜在风险与解决策略

zadig安装驱动潜在风险与解决策略 ✨没事不要闲着乱打驱动&#xff0c;能正常使用的情况下&#xff0c;不要轻易或随意去乱打驱动&#xff0c;可能会导致新的驱动对已有的设备不兼容的问题。✨&#x1f530;特别说明&#xff1a;本文介绍的方法&#xff0c;并不能包治百病&…

Android学习之路(2) 文本设置

Android学习之路(1) 文本 一、设置文本内容 设置文本内容的两种方式&#xff1a; 一种是在XML文件中通过属性android:text设置文本代码如下 <TextViewandroid:id"id/tv_hello"android:layout_width"wrap_content"android:layout_height"wrap_c…

20套面向对象程序设计选题Java Swing(含教程)持续更新

20套面向对象程序设计选题&#xff0c;适合Java课程设计&#xff0c;可用MySQL数据库&#xff0c;也可以不使用数据库&#xff0c;使用Java集合存储数据。 持续更新&#xff0c;建议收藏 点击获取代码 0. JavaSwing管理系统万能模板 视频教程&#xff1a; 【课程设计】2小时…

报错Uncaught (in promise) Error: Manifest request to...

在使用nuxt框架时&#xff0c;出现如下报错&#xff1a; 解决方案&#xff1a; 不要打开两个以上的开发者工具更换nuxt的端口号 参考资料&#xff1a;https://github.com/nuxt/nuxt.js/issues/6202

将tp5项目、fastadmin项目部署到服务器宝塔面板

目录 一、将你的fastadmin或者tp5项目文件夹上传至你的服务器域名根目录下 二、修改你的网站目录指向&#xff0c;指向public目录&#xff0c;点击保存&#xff0c;并取消勾选防跨站攻击。 三、配置伪静态 四、fastadmin框架上传至服务器后如果想要访问后台可以进行重定向&am…

低成本NFC端口静电保护方案图及ESD二极管选型指南

Near Field Communication&#xff0c;简称&#xff1a;NFC&#xff0c;中文名称&#xff1a;近场通信&#xff0c;是一种短距离高频的无线电技术&#xff0c;能够实现近距离无线通讯和数据交换&#xff0c;是由非接触式射频识别&#xff08;RFID&#xff09;及互连互通技术整合…

‘float‘ object has no attribute ‘decode‘

最近在jieba分词的时候报错&#xff1a;float object has no attribute decode 我的源代码如下 data pd.read_csv(../data/answers_covid.csv,encodingutf-8) 解决方法:.astype(str) data pd.read_csv(../data/answers_covid.csv,encodingutf-8).astype(str)

2023/8/9总结

首页数据分页&#xff1a; 消息列表的完善&#xff1a; 文章水印&#xff1a;

【Vue3】keep-alive 缓存组件

当在 Vue.js 中使用 <keep-alive> 组件时&#xff0c;它将会缓存动态组件&#xff0c;而不是每次渲染都销毁和重新创建它们。这对于需要在组件间快速切换并且保持组件状态的情况非常有用。 <keep-alive> 只能包含&#xff08;或者说只能渲染&#xff09;一个子组件…

使用chatGPT生成提示词,在文心一言生成装修概念图

介绍 家是情感的港湾&#xff0c;而家居装修则是将情感融入空间的艺术。如何在有限的空间里展现个性与美感&#xff0c;成为了现代人关注的焦点。而今&#xff0c;随着人工智能的发展&#xff0c;我们发现了一个新的创意助手——ChatGPT&#xff0c;它不仅为我们带来了更多可能…

C++:类与对象(上):类结构体概念、访问限定符、对象实例化、类对象模型

目录 一、认识面向过程和面向对象 1.1 面向过程 1.2 面向对象 二、初步了解类与结构体 2.1 结构体&#xff08;struct&#xff09; 2.2 类&#xff08;class&#xff09; 三、类 3.1 语法 3.2 定义方式 3.2.1 不分文件编写 3.2.2 分文件编写&#xff08;常用&#…

SpringBoot案例-部门管理-查询

查看页面原型&#xff0c;明确需求需求 页面原型 需求分析 阅读接口文档 接口文档链接如下&#xff1a; https://onedrive.live.com/?cidC62793E731F0C1BE&idC62793E731F0C1BE%2148 思路分析 用户发送请求&#xff0c;交由对应的Controller类进行处理&#xff0c;Con…

初识Java集合框架

前言 在大多数情况下&#xff0c;你常常会看到《C数据结构》类似的书籍&#xff0c;很多人可能会认为数据结构是一门依赖语言的学科&#xff0c;有了语言才可能有数据结构&#xff0c;其实这里需要纠正的是&#xff0c; 数据结构(Data Structure)是计算机存储、组织数据的方式…

二、web核心防御机制(上)

文章目录 一、Web应用程序与风险1.1 Web应用程序的发展历程1.2 Web应用程序安全 二、核心防御机制2.1 处理用户访问2.2 处理用户输入2.2.1输入的多样性2.2.2 输入处理方法2.2.3 边界确认2.2.4 多步确认与规范化 一、Web应用程序与风险 1.1 Web应用程序的发展历程 在因特网发展…

Xception 算法详解

📮 本次重点(模型轻量化): ● Inception设计理念 ● 点卷积 ● 深度可分离卷积 ● Bottleneck结构 注:Xception算法整体结构是其次,主要是了解以上四个结构。 今天详解Xception算法,由于Xception模型在极大的减少了网络参数量和计算复杂度的同时,可以保持卓越的性能表…

优思学院|精益生产3大特征:拉动式生产、消除浪费、自働化

在现代制造业中&#xff0c;精益生产作为一种高效的生产方式&#xff0c;受到了广泛的关注和应用。它的核心理念是通过最大程度地减少浪费和提高效率来实现生产过程的优化。经典的精益体系从中提炼出了精益生产的三大特征&#xff0c;分别是拉动式生产、消除浪费以及自働化。让…