matlab中的双层数值积分
🏆本文收录于《CSDN问答解惑-专业版》专栏，主要记录项目实战过程中的Bug之前因后果及提供真实有效的解决方案，希望能够助你一臂之力，帮你早日登顶实现财富自由🚀；同时，欢迎大家关注&&收藏&&订阅！持续更新中，up！up！up！！
问题描述
matlab中的双层数值积分.在matlab中复现一篇文献的公式和结果，但是跑出来的图像和文献相比向高温区偏移了120℃左右，帮我看看哪儿有问题。
涉及的公式：da/dt=\int_{0}{\infty}\frac{A}{\beta}exp\left[-\frac{E}{RT}-\frac{A}{\beta}\int_{T_0}{T}exp\left(-\frac{E}{RT}\right)dT\right]f\left(E\right)dE（latex格式）
其中f（E）是高斯分布，内层温度积分由于没有解析解，所以用了其他论文的近似式
文献：An overview of distributed activation energy model and its application
in the pyrolysis of lignocellulosic biomass.式（17）（18），复现图1
%% 传统分布式活化能计算
clc;clear all;
% 常量设置
R = 8.314;
beta = 10; % 升温速率
b = 10;% 温度间隔A1 = 1.67e13; % 频率因子
A2 = 1e15;
A3 = 1e16;
E01 = 228000; % 平均活化能% 积分上下限设置
Ta = 500;
Tb = 1000;
j = (Tb-Ta)/b;
x = linspace(Ta,Tb,j+1);figure1 = figure('Name','高斯拟合da_dt');
% figure2 = figure('Name','高斯拟合a_t');tao = 0.0005; % 活化能分布的标准差相对于E0的百分比
sigma1 = tao*E01; % 活化能分布的标准差，J/mol
for tao = 0.005:0.01:0.1sigma1 = tao*E01;f1 = @(E) normpdf(E,E01,sigma1);% E_scope =linspace(0,E01+30*sigma1,1000); % 温度积分近似式dI1 = @(E,T) (((R*T*T)./E).*(((0.99962.*E+0.60642*R*T)/(E+2.56879*R*T)).*exp(-E./(R*T))));%-((R*Ta*Ta)./E).*(((0.99962.*E+0.60642*R*Ta)./(E+2.56879*R*Ta)).*exp(-E./(R*Ta))));da_scope1 = zeros(j,1);a_scope1 = zeros(j,1);% 整体积分结果da_dt1 = @(E,T) (A1/beta)*exp((-E./(R*T))-((A1/beta))*dI1(E,T)).*f1(E);
% a_t1 = @(E,T) (A1/beta).*exp(-dI1(E,T)).*f1(E);k = 0;for n = Ta:b:Tbk = k+1;
% a_scope1(k) = 1-integral(@(E) a_t1(E,n),0,E01+30*sigma1);da_scope1(k) = integral(@(E) da_dt1(E,n),0,E01+30*sigma1);endfigure1;plot(x,da_scope1);hold on;
% figure2;
% plot(x,a_scope1);
% hold on;
end
如上问题有来自我自身项目开发，有的收集网站，有的来自读者，如有侵权，立马删除。
解决方案
如下是上述问题的解决方案，仅供参考：
在您提供的 MATLAB 代码中，您尝试复现文献中的双层数值积分问题，但结果与文献相比有偏差。以下是一些可能导致问题的原因和相应的建议：
-
积分上下限：在您的代码中，内层积分的上下限设置可能不正确。通常，内层积分的上限应该是无穷大，而不是
E01+30*sigma1。您可能需要调整这个上限，或者使用 MATLAB 的integral函数来处理无限区间的积分。 -
温度积分近似式：您使用了文献中的近似式
dI1来计算内层积分。请确保这个近似式是正确的，并且适用于您的问题。如果这个近似式不正确或不适用于您的数据，它可能会导致结果偏差。 -
高斯分布函数：确保
normpdf函数的参数设置正确，特别是均值E01和标准差sigma1。 -
积分函数：在定义
da_dt1函数时，您使用了integral函数来计算积分。请确保integral函数的调用是正确的，并且积分的上下限是适当的。 -
温度范围：在循环中，您使用了
for n = Ta:b:Tb来遍历温度范围。这个循环可能没有正确地覆盖整个温度范围。您可能需要调整这个循环的起始点和步长。 -
绘图：在您的代码中，您使用了
plot函数来绘制结果。请确保您的x向量包含了正确的温度范围，并且da_scope1向量包含了在该温度范围内计算的积分结果。 -
文献中的公式：确保您的 MATLAB 代码正确地实现了文献中的公式。如果可能，请与文献作者联系，获取更多的实现细节或验证代码。
-
数值方法的准确性：数值积分方法（如
integral函数）可能会受到数值误差的影响。您可以尝试增加积分的精度或使用不同的数值积分方法来检查结果的稳定性。 -
代码逻辑：检查代码中的逻辑是否正确，特别是涉及到循环和条件判断的部分。
-
环境设置：确保您的 MATLAB 环境设置正确，没有其他变量或函数影响结果。
如果上述建议仍然无法解决问题，您可能需要进一步检查文献中的公式和实现细节，或者与领域专家讨论以获得帮助。此外，您也可以尝试使用不同的数值方法或软件来验证您的结果。
希望如上措施及解决方案能够帮到有需要的你。
PS：如若遇到采纳如下方案还是未解决的同学，希望不要抱怨&&急躁，毕竟影响因素众多，我写出来也是希望能够尽最大努力帮助到同类似问题的小伙伴，即把你未解决或者产生新Bug黏贴在评论区，我们大家一起来努力，一起帮你看看，可以不咯。
若有对当前Bug有与如下提供的方法不一致，有个不情之请，希望你能把你的新思路或新方法分享到评论区，一起学习，目的就是帮助更多所需要的同学，正所谓「赠人玫瑰，手留余香」。
☀️写在最后
ok，以上就是我这期的Bug修复内容啦，如果还想查找更多解决方案，你可以看看我专门收集Bug及提供解决方案的专栏《CSDN问答解惑-专业版》，都是实战中碰到的Bug，希望对你有所帮助。到此，咱们下期拜拜。
码字不易，如果这篇文章对你有所帮助，帮忙给 bug菌 来个一键三连(关注、点赞、收藏) ，您的支持就是我坚持写作分享知识点传播技术的最大动力。
同时也推荐大家关注我的硬核公众号:「猿圈奇妙屋」 ；以第一手学习bug菌的首发干货，不仅能学习更多技术硬货，还可白嫖最新BAT大厂面试真题、4000G Pdf技术书籍、万份简历/PPT模板、技术文章Markdown文档等海量资料，你想要的我都有！
📣关于我
我是bug菌，CSDN | 掘金 | InfoQ | 51CTO | 华为云 | 阿里云 | 腾讯云 等社区博客专家，C站博客之星Top30，华为云2023年度十佳博主，掘金多年度人气作者Top40，掘金等各大社区平台签约作者，51CTO年度博主Top12，掘金/InfoQ/51CTO等社区优质创作者；全网粉丝合计 30w+；硬核微信公众号「猿圈奇妙屋」，欢迎你的加入！免费白嫖最新BAT互联网公司面试真题、4000G PDF电子书籍、简历模板等海量资料，你想要的我都有，关键是你不来拿哇。
