Question

MATLAB编程紧急求助

function y=K(k)
%program to compute the complete eliptic integrals K(k)
%Using the arith-geom mean method
[m,n]= size(k);
if m+ n==2
k=[k k];
end
if m==1
k= k';
end
a(:, 1)=ones(k);
b(:, 1)=sqrt(ones(k)-k.*k);
c(:, 1)=k;
i=1;
while max(c(: ,1))>eps
a(:, i+ 1)=.5*(a(:,i)+b(:,i));
b(:, i+ 1)=sqrt(a(:,i).*b(:,i));
c(:, i+ 1)=.5*(a(:,i)-b(:,i));
i=i+1;
end
t=pi*ones(k)./(2*a(:,i));
if m+n==2
y=t(1,:);
elseif m==1
y= t';
else
y= t;
end

运行输入 K（[1 2 3 4 ]）

提示：

??? Error using ==> ones
Size vector must be a row vector with integer elements.

Error in ==> K at 11
a(:, 1)=ones(k);

怎么回事呢

Answer 1

哎，老了，matlab生疏了。
哥们啊，我一个心血来潮给你认真了一个下午才弄出来的，要求追加分数！
程序代码如下：
beta=0.1;
w0=1;
delta_t=25/1000;
t(1)=0;
for i=1:1000
t(i+1)=delta_t*i;
end

threshold=0.9;
x(1)=100; %为了求解方便，最后再更改比例尺
temp_v1=0;
n=1; %x的角标
for i=1:1000
for j=-100:0.001:100
temp_x=j;
temp_v2=(temp_x-x(n))/delta_t;
temp_a=(temp_v2-temp_v1)/delta_t;
if abs(temp_a+2*beta*temp_v2+w0^2*temp_x)<threshold
n=n+1;
x(n)=temp_x;
temp_v1=temp_v2;
break;
end
end
end
for i=1:1001 %回归真正的比例
x(i)=x(i)/100;
end
plot(t,x);
xlabel('Time');
ylabel('X-distance');
title('弹簧振子阻尼振动曲线');

直接运行便可以得到弹簧振子阻尼振动曲线。

Answer 2

问题出在这里 ：
a(:, 1)=ones(k);
b(:, 1)=sqrt(ones(k)-k.*k);
c(:, 1)=k;
ones（k）生成的是一个全1的矩阵，a(:, 1)是一列数据 维数不匹配