|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
%product.m%
function y=product(a0)
a=a0;
if (a<21)
y=factorial(a);
else
y=1;
for i=1:a
y=y*i
end
end
%sun.m%
function y=sun(k0,t0,w10,w20,x)
k=k0;t=t0;w1=w10;w2=w20;
y=0;
for p=0:k
sum1=0;n=p;
while( n<=2*t)
if (n~=t)
s=(sin(w2*(n-t))-sin(w1*(n-t)))/(pi*(n-t));
sum1=sum1+s*product(k+n-p)/(product(n-p)*product(k))*x.^(k+n-2*p);
else
s=(w2-w1)/pi;
sum1=sum1+s*product(k+n-p)/(product(n-p)*product(k))*x.^(k+n-2*p);
end
n=n+1;
end
z=(-1)^(p+k)*product(k)/(product(k-p)*product(p));
y=y+z*sum1;
end
%fun4.m%
function y=fun4(x)
y=sun(20,16,3*pi/4,pi,x)
%newton2.m%
function y=newton2(x01,x02,m,eps)
x(1)=x01;x(2)=x02;
c=1;
i=2;
while (abs(c)>eps*x(i))
x(i+1)=x(i)-fun4(x(i))*(x(i)-x(i-1))/(fun4(x(i))-fun4(x(i-1)));
c=x(i+1)-x(i);
i=i+1;
if (i>m)error('m is full');
return;
end
end
y=x(i);i
从图中可以看出其x在(-1,1)之间存在很多值,但是当我设定初值为负值时,结果经常是NaN,有时是被零除,有时候没有任何说明和警告。请高手指教,谢谢! |
|