查看: 93|回复: 1|关注: 0

[已答复] 请问如图所示公式应该用哪种数值积分

[复制链接]

新手

5 麦片

财富积分


050


2

主题

2

帖子

0

最佳答案
发表于 2019-11-6 15:01:51 | 显示全部楼层 |阅读模式
本帖最后由 demonseasons 于 2019-11-6 15:11 编辑

如图所示,公式中包含累加和积分,累加我使用for循环,但数值积分部分使用哪个积分搞不清楚,目前已用过trapz和quad、quadl,得到结果均不对
公式所对应代码如下:
  1. for m1=1:Com
  2.     for j=1:size(tD,2)  %时间循环
  3.         f=log(2)/tD(j);
  4.         PwD(m1,j)=0;
  5.         for i=1:N
  6.             z=i*f;
  7.               beta=n*pi;
  8.               kesai0=sqrt(z*hD^2*LD^2+beta.*LD^2);
  9.               kesai1=besselk(0,sqrt((xD-arph).^2).*kesai0);   %被积分函数
  10.               int5=besselk(0,sqrt((xD-arph).^2).*sqrt(z*hD^2*LD^2));
  11.               int6=@(arph)besselk(0,sqrt((xD-arph).^2).*sqrt(z*hD^2*LD^2));
  12.               int3=quadgk(int6,-1,1);
  13.               int1=@(arph)besselk(0,sqrt((xD-arph).^2).*kesai0).*cos(beta*zrD).*cos(beta*zwD);
  14.               int2=0;
  15.               for n=1:1:1000       %累加函数
  16.               int2=int2+quadgk(int1,-1,1);
  17.               end
  18. pD(m1,j)=0.5./z.*(int3+2.*int2);
  19.                 PwD(m1,j)=1/(z*(z+1/(z*pD(m1,j)+hD*S)));
  20.                 PwD(m1,j)=PwD(m1,j)+V(i)*pD(m1,j);
  21.         end
  22.             
  23.             PwD(m1,j)=PwD(m1,j)*f;
  24.     end
  25. end   
复制代码


m文件如下,烦请各位老师指正



公式

公式

horizontal.m

2.71 KB, 下载次数: 0

新手

19 麦片

财富积分


050


8

主题

39

帖子

2

最佳答案
发表于 2019-11-6 15:29:48 | 显示全部楼层
涉及各种积分变换  ,,这个时候matlab的积分会偶尔不准,,这是正常现象。说明积分过程中有不可导或导不连续的点,用哪个都不行
建议如果是拉普拉斯变换的话  不要用定义法,,直接用 laplace(f1,t,s)试试
逆变换是ilaplace
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表