查看: 292|回复: 5|关注: 0

[已解决] 用梯形法计算一个二重积分,运行结果是NaN,能看看哪出问题了吗?

[复制链接]

新手

7 麦片

财富积分


050


2

主题

6

帖子

0

最佳答案
本帖最后由 Gal. 于 2020-3-29 21:41 编辑

function s=dou_integral(func,a,b,c,d)
n=100;
err=1e-4;
hx=(b-a)/n;
hy=(d-c)/n;
s=0;
for k=1:n-1   
     x=a+hx*k;
    for t=1:n-1
        y=c+hy*t;
        s=s+feval(func,x,y);
    end
end
s2=0;
for k=1:n-1
    y=c+hy*k;
s2=s2+feval(func,a,y)+feval(func,b,y);
end
s3=0;
for k=1:n-1
    x=c+hx*k;
s3=s3+feval(func,x,c)+feval(func,x,d);
end
s1=hx*hy*((feval(func,d,b)+feval(func,b,c)+feval(func,c,a)+feval(func,d,a))/4+s+(s2+s3)/2);

n=n*2;
hx=(b-a)/n;
hy=(d-c)/n;
s=0;
for k=1:n-1   
     x=a+hx*k;
    for t=1:n-1
        y=c+hy*t;
        s=s+feval(func,x,y);
    end
end
s2=0;
for k=1:n-1
    y=c+hy*k;
s2=s2+feval(func,a,y)+feval(func,b,y);
end
s3=0;
for k=1:n-1
    x=c+hx*k;
s3=s3+feval(func,x,c)+feval(func,x,d);
end
s=hx*hy*((feval(func,d,b)+feval(func,b,c)+feval(func,c,a)+feval(func,d,a))/4+s+(s2+s3)/2);
%循环计算达到精度
while abs(s-s1)>err

    s1=s;
    n=n*2;
    hx=(b-a)/n;
hy=(d-c)/n;
s=0;
for k=1:n-1   
     x=a+hx*k;
    for t=1:n-1
        y=c+hy*t;
        s=s+feval(func,x,y);
    end
end
s2=0;
for k=1:n-1
    y=c+hy*k;
s2=s2+feval(func,a,y)+feval(func,b,y);
end
s3=0;
for k=1:n-1
    x=c+hx*k;
s3=s3+feval(func,x,c)+feval(func,x,d);
end
s=hx*hy*((feval(func,d,b)+feval(func,b,c)+feval(func,c,a)+feval(func,d,a))/4+s+(s2+s3)/2);
end
s  %#ok<NOPRT>

这是运行结果

这是运行结果
回复主题 已获打赏: 0 积分

举报

新手

7 麦片

财富积分


050


2

主题

6

帖子

0

最佳答案
 楼主| 发表于 2020-3-29 21:35:50 | 显示全部楼层
补充一下func
func=@(x,y) sin(x+y)./(x+y)
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

3

主题

1261

帖子

272

最佳答案
  • 关注者: 58
发表于 2020-3-29 22:16:25 | 显示全部楼层 |此回复为最佳答案
运行这句代码s=s+feval(func,x,y);时,x+y有时等于0,即分母为零,导致s不是数,即NaN
回复此楼 已获打赏: 0 积分

举报

新手

7 麦片

财富积分


050


2

主题

6

帖子

0

最佳答案
 楼主| 发表于 2020-3-29 22:22:25 | 显示全部楼层
20141303 发表于 2020-3-29 22:16
运行这句代码s=s+feval(func,x,y);时,x+y有时等于0,即分母为零,导致s不是数,即NaN ...

你好,请问怎样解决呢
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

3

主题

1261

帖子

272

最佳答案
  • 关注者: 58
发表于 2020-3-30 09:04:59 | 显示全部楼层
建议针对积分区间做个预判,如果超出定义域,就输出警告,避免出现这种情况,或者用个比较小的数代替0,当然分寸不太好把握,仅供参考
回复此楼 已获打赏: 0 积分

举报

新手

7 麦片

财富积分


050


2

主题

6

帖子

0

最佳答案
 楼主| 发表于 2020-3-30 17:08:31 | 显示全部楼层
20141303 发表于 2020-3-30 09:04
建议针对积分区间做个预判,如果超出定义域,就输出警告,避免出现这种情况,或者用个比较小的数代替0,当 ...

谢谢:loveliness:
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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