请用MATLAB画含有二重积分和累加求和的随时间t变化的函数图像，不知如何把二重积分加到循环里，请解答一下

Question

想要画一个关于时间t变化的图像，由两部分时间t变化的函数组成，其中一个式子P1包含了二重积分求解，我只写了循环求和，对于有变量的二重积分求解不知如何加进去，请解答一下，谢谢！

公式如下：

代码如下：

close all;
clear;
clc;

D = 0.00125;
L = 1.5;
muw = 1.2;
sgmw = 0.2;
mub = 8.4823e-9;
sgmb = 6.0016e-10;
fai = 0;
muy = 1e-4;
sgmy = 2e-5;
lmd = 2.5e-5;

tx = 0:1:5000;
p1 = 9.9931015E-01 + 6.8984636E-04*exp(-1.5510700E+00*tx);
p2 = 9.9993146E-01 + 6.8538349E-05*exp(-1.5101035E+00*tx);
p3 = 9.9993400E-01 + 6.5995644E-05*exp(-1.5501023E+00*tx);
p4 = 9.9993633E-01 + 6.3669415E-05*exp(-1.4700936E+00*tx);

n_max = 100; % 设置累加求和的最大次数

% 初始化结果数组
r2 = zeros(size(tx));
r3 = zeros(size(tx));
r4 = zeros(size(tx));
G5 = zeros(size(tx));

% 第一部分
m1 = normcdf(D, muw, sgmw);

for n = 0:n_max
    % 第二部分,二重积分
   
    
    % 第三部分
    m3 = exp(-lmd .* tx) .* ((lmd .* tx) .^ n) ./ factorial(n);
    
    % 累加求和
    r1 = m1 .^ n .* m2 .* m3;
    r2 = r2 + r1;
end

% 计算1-r2
r3 = 1 - r2;
G5 = (1-(1-p1).*(1-p2)).* p3.*p4;

r4 = r3.*G5;
plot(tx, r4, 'k', 'LineWidth', 1.5);
set(gca, 'FontSize', 14, 'Fontname', 'Times New Roman')
xlabel("Time/Day", 'fontsize', 14, 'FontName', 'Times New Roman', 'fontweight', 'bold');
ylabel("R", 'fontsize', 14, 'FontName', 'Times New Roman', 'fontweight', 'bold', 'Fontangle', 'italic');

Answer 1

该回答引用讯飞星火及结果代码方法已验证可行
根据您提供的公式和代码，我理解您想要将二重积分求解添加到循环求和中。为了实现这个目标，您可以使用MATLAB中的integral函数来计算二重积分。以下是修改后的代码：

close all;
clear;
clc;
 
D = 0.00125;
L = 1.5;
muw = 1.2;
sgmw = 0.2;
mub = 8.4823e-9;
sgmb = 6.0016e-10;
fai = 0;
muy = 1e-4;
sgmy = 2e-5;
lmd = 2.5e-5;
 
tx = 0:1:5000;
p1 = 9.9931015E-01 + 6.8984636E-04*exp(-1.5510700E+00*tx);
p2 = 9.9993146E-01 + 6.8538349E-05*exp(-1.5101035E+00*tx);
p3 = 9.9993400E-01 + 6.5995644E-05*exp(-1.5501023E+00*tx);
p4 = 9.9993633E-01 + 6.3669415E-05*exp(-1.4700936E+00*tx);
 
n_max = 100; % 设置累加求和的最大次数
 
% 初始化结果数组
r2 = zeros(size(tx));
r3 = zeros(size(tx));
r4 = zeros(size(tx));
G5 = zeros(size(tx));
 
% 第一部分
m1 = normcdf(D, muw, sgmw);
 
for n = 0:n_max
    % 第二部分，二重积分
    integral_result = integral(@(x) (normpdf(x, mub, sgmb)), -Inf, Inf);
    m2 = integral_result(1);
    
    % 第三部分
    m3 = exp(-lmd .* tx) .* ((lmd .* tx) .^ n) ./ factorial(n);
    
    % 累加求和
    r1 = m1 .^ n .* m2 .* m3;
    r2 = r2 + r1;
end
 
% 计算1-r2
r3 = 1 - r2;
G5 = (1-(1-p1).*(1-p2)).* p3.*p4;
 
r4 = r3.*G5;
plot(tx, r4, 'k', 'LineWidth', 1.5);
set(gca, 'FontSize', 14, 'Fontname', 'Times New Roman')
xlabel("Time/Day", 'fontsize', 14, 'FontName', 'Times New Roman', 'fontweight', 'bold');
ylabel("R", 'fontsize', 14, 'FontName', 'Times New Roman', 'fontweight', 'bold', 'Fontangle', 'italic');

这段代码将在循环求和中添加二重积分求解的部分。希望这对您有所帮助！