matlab的参数积分编程

就是h的变化是从0到400，需要求的东西是，x和y在这段过程的所描绘出来的曲线。这个要怎么编程啊！

举报该问题

推荐答案推荐于2016-05-22

前几天写了这个问题的程序，因时间问题未能及时回答。现看到楼上已经回答，但代码以收费附件的方式提供，不便于交流，我谈谈自己的做法吧。

注意：以下列出几种方法的代码，是各自独立的，可以分别直接复制到命令窗口里面运行，或者保存成M文件之后运行。

方法1：符号积分

这种方法最为简单直接，很容易看懂。使用符号运算求积分，对于无法求出解析解的情况，有时候可以计算得到高精度数值解。

参考代码如下：

syms h
phi = atan(h/10.8);
alpha = int(sin(phi),0,h)/210;
T = linspace(0,400,20);
X = T*0;
Y = X;
for i=1:length(T)
    t = T(i);
    x = int(37.5*cos(phi+alpha),0,t);
    X(i) = double(x);
    y = int(37.5*sin(phi+alpha),0,t);
    Y(i) = double(y);
end
plot(X,Y)
axis equal
xlabel x; ylabel y

在有些版本中会提示“无法求出显式解”（Explicit integral could not be found）的警告信息，如果不希望看到警告，可以在代码前加一句

warning off symbolic:sym:int:warnmsg1

方法2：数值积分

这种方法通用性更强，对于一些很复杂的积分，使用符号运算可能会失败，而用数值方法一般不会有问题。

参考代码如下：

quadfun = @quadgk;
phi = @(h)atan(h/10.8);
alpha = @(H)arrayfun(@(h)quadfun(@(x)sin(phi(x)),0,h)/210,H);
H = linspace(0,2000,50);
X = arrayfun(@(h)quadfun(@(x)37.5*cos(phi(x)+alpha(x)),0,h),H);
Y = arrayfun(@(h)quadfun(@(x)37.5*sin(phi(x)+alpha(x)),0,h),H);
plot(X,Y)
axis equal
xlabel x; ylabel y

说明几点：

1、这里涉及到的知识主要包括匿名函数、arrayfun和数值积分函数，如果看不明白，可以针对这几个方面查阅相关资料。

2、为了更清楚地了解x-y之间的关系，把h的取值范围适当放大了一些（0~2000）。

3、关于积分函数的选择：个人认为，使用quadgk进行积分是一个比较好的选择（该函数需2007b或更高版本）。如果使用2012a之后的版本，也可以用integral函数。如果是早期版本的MATLAB，可以用quadl之类的函数，但计算效率会低很多，而且可能产生关于最小步长的警告（可以用warning off MATLAB:quadl:MinStepSize关掉）。

关于补充问题

很明显，新图片里面的问题与原问题不同，其它系数是否一致姑且不论，最大的差别是，第一个式子原来是关于h的反正切函数，而现在是关于sin(h)的反正切。这个差别可能会导致方法1（使用符号运算）失效，而使用方法2（数值积分）求解的时间也远远超过原问题（就我这里测试的情况看，所需时间大约是40倍）。

有个系数R没有给出，我按照原问题对应的系数设置 R = 37.5*r0*i1 = 112.5000，参考代码如下（h取值范围0~2000，可自行修改）：

r3 = 15; C = 33;
r0 = 10; i1 = 0.3;
i2 = 2; l = 90;
m = 190;
R = 37.5*r0*i1;
quadfun = @quadgk;
phi = @(h) atan( r3*sin(h/(r0*i1*i2)) / C );
alpha = @(H) arrayfun(@(h)quadfun(@(x)sin(phi(x)),0,h)/sqrt(l^2+m^2),H);
H = linspace(0,2000,50);
%x = quadl(@(x)37.5*cos(phi(x)+alpha(x)),0,1);
X = arrayfun(@(h)quadfun(@(x)R/(r0*i1)*cos(phi(x)+alpha(x)),0,h),H);
Y = arrayfun(@(h)quadfun(@(x)R/(r0*i1)*sin(phi(x)+alpha(x)),0,h),H);
plot(X,Y)