求超越方程的根
X=-0.1:0.001:0.1;a=130;
b=200;
J00=besselj(n,130*X);
J01=besselj(n+1,130*X);
J10=besselj(n,200*X);
J11=besselj(n+1,200*X);
Y00=bessely(n,130*X);
Y01=bessely(n+1,130*X);
Y10=bessely(n,200*X);
Y11=bessely(n+1,200*X);
K=(n/a-1/5000)*J00-X.*J01;
L=(n/b+1/5000)*Y10-X.*Y11;
V=(n/b+1/5000)*J10-X.*J11;
W=(n/a-1/5000)*Y00-X.*Y01;
y=K.*L-V.*W;
plot(X,y)
hold on
plot([-0.15 0.15],,'k');
f=@(X)K.*L-V.*W;
X=[-0.15 -0.1];
arrayfun(@(X)fzero(f,X),X);
Warning: Imaginary parts of complex X and/or Y arguments ignored
??? Operands to the || and && operators must be convertible to logical scalar values.
Error in ==> fzero at 323
elseif ~isfinite(fx) || ~isreal(fx)
Error in ==> @(X)fzero(f,X) 你定义的 f 并非是一个含有未知数 X 的函数句柄，而是确定的数值，因为 K.*L-V.*W 被计算为确定的数了。按照你目前的写法，可以将 J00 ~ J11, Y00~Y11，K, L, V, W 全部定义为关于X的函数，后面用
f = @(X)K(X).*L(X)-V(X).*W(X);
另外，你的 n 是多少？ 补充一下，你的方程并不是很复杂，其实没有必要定义这么多中间变量 J00 ~ J11, Y00~Y11。如果你是担心定义 f 的时候表达式过长， 只需分别定义 K, L, V, W 函数，然后直接定义 f 谢谢！我试试！ 我试了一下，还是不行。
n=0;
X=-0.1:0.001:0.1;
a=130;
b=200;
y=((n/a-1/5000)* besselj(n,130*X)-X.* besselj(n+1,130*X)).* ((n/b+1/5000)* bessely(n,200*X)-X.* bessely(n+1,200*X))- ((n/b+1/5000)* besselj(n,200*X)-X.* besselj(n+1,200*X)).* ((n/a-1/5000)* bessely(n,130*X)-X.* bessely(n+1,130*X));
plot(X,y)
hold on
plot([-0.15 0.15],,'k');
f=@(X) ((n/a-1/5000)* besselj(n,130*X)-X.* besselj(n+1,130*X)).* ((n/b+1/5000)* bessely(n,200*X)-X.* bessely(n+1,200*X))- ((n/b+1/5000)* besselj(n,200*X)-X.* besselj(n+1,200*X)).* ((n/a-1/5000)* bessely(n,130*X)-X.* bessely(n+1,130*X));
X0=-0.09;
arrayfun(@(X)fzero(f,X),X0);
Warning: Imaginary parts of complex X and/or Y arguments ignored
??? Error using ==> fzero at 324
Function value at starting guess must be finite and real.
Error in ==> @(X)fzero(f,X)
wuhaimiao 发表于 2014-1-26 18:56
我试了一下，还是不行。
n=0;
X=-0.1:0.001:0.1;
n=0;
X=-0.1:0.001:0.1;
a=130;
b=200;
f=@(X) ((n/a-1/5000)* besselj(n,130*X)-X.* besselj(n+1,130*X)).* ((n/b+1/5000)* bessely(n,200*X)-X.* bessely(n+1,200*X))- ((n/b+1/5000)* besselj(n,200*X)-X.* besselj(n+1,200*X)).* ((n/a-1/5000)* bessely(n,130*X)-X.* bessely(n+1,130*X));
fplot(f,[-0.1,0.1])
hold on
plot([-0.15 0.15],,'k');
X0=;
sol = arrayfun(@(X)fzero(f,X),X0)
sol = ;
plot(sol, f(sol), 'bo')
万分感谢！ 万分感谢！
页:
[1]