查看: 297|回复: 2|关注: 0

[已解决] matlab用trapzq求解数值积分出错

[复制链接]

新手

10 麦片

财富积分


050


5

主题

12

帖子

0

最佳答案
微信图片_20191109190611.jpg 其中F(k)和E(k)分别为第一类和第二类完全椭圆积分
求M21,已知
miu=(4*pi)*1e-7;
d=30.48;
Rs=15.24;
Rp=15.24;
c=0
下面是我的程序:
clc;
clear;
%h=0:pi;
miu=(4*pi)*1e-7;
d=30.48;
Rs=15.24;
Rp=15.24;
%N=2*miu*sqrt(Rs*Rp)/pi;
c=0;
a=Rs/Rp;
b=c/Rp;
syms h
V=sqrt(1+d^2/Rs^2-2*d*cos(h)/Rs);
k=sqrt(4*a*V/((1+a*V)^2+b^2));
[K,E] = ellipke(k^2);
P=@(k)((1-k^2/2)*K-E);
f=P(k)*(1-d*cos(h)/Rs)/(k*(sqrt(V^3)));
%F=int(f,0,pi);
h=0:pi/100:pi;
% g=length(f);
F=trapz(h,f,2)
ans=vpa(F)

报错:
搜狗截图19年11月09日1912_1.jpg

论坛优秀回答者

0

主题

527

帖子

100

最佳答案
  • 关注者: 20
发表于 2019-11-10 20:40:53 | 显示全部楼层 |此回复为最佳答案
不应该使用syms h,该位置应把后面的h=0:pi/100:pi;换过来,我简单改了一下,具体需要你自己改
clc;
clear;
%h=0:pi;
miu=(4*pi)*1e-7;
d=30.48;
Rs=15.24;
Rp=15.24;
%N=2*miu*sqrt(Rs*Rp)/pi;
c=0;
a=Rs/Rp;
b=c/Rp;
h=0:pi/100:pi;
V=sqrt(1+d^2/Rs^2-2*d*cos(h)/Rs);
k=sqrt(4*a.*V./((1+a.*V).^2+b^2));
[K,E] = ellipke(k.^2);
P=@(k)((1-k.^2./2).*K-E);
f=P(k).*(1-d.*cos(h)./Rs)./(k.*(sqrt(V.^3)));
f(1)=0;
%F=int(f,0,pi);
h=0:pi/100:pi;
% g=length(f);
F=trapz(h,f,2);
ans=vpa(F)

结果如图
QQ图片20191110203930.png

新手

10 麦片

财富积分


050


5

主题

12

帖子

0

最佳答案
 楼主| 发表于 2019-11-14 21:33:11 | 显示全部楼层
大哥  太棒了  谢谢宁!帮我解决了一个大问题
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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