求超越方程的根

wuhaimiao

  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],[0 0],'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)

winner245

你定义的 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 是多少？

winner245

补充一下，你的方程并不是很复杂，其实没有必要定义这么多中间变量 J00 ~ J11, Y00~Y11。如果你是担心定义 f 的时候表达式过长， 只需分别定义 K, L, V, W 函数，然后直接定义 f

wuhaimiao

谢谢！我试试！

wuhaimiao

我试了一下，还是不行。
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],[0 0],'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)

winner245

wuhaimiao 发表于 2014-1-26 18:56
我试了一下，还是不行。
n=0;
X=-0.1:0.001:0.1;

1. n=0;
   
2. X=-0.1:0.001:0.1;
   
3. a=130;
   
4. b=200;
   
5. 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));
   
6. fplot(f,[-0.1,0.1])  
   
7. hold on  
   
8. plot([-0.15 0.15],[0 0],'k');
   
9. X0=[0.09 0.05 0.01];
   
10. sol = arrayfun(@(X)fzero(f,X),X0)
    
11. sol = [sol,-sol];
    
12. plot(sol, f(sol), 'bo')

复制代码

wuhaimiao

万分感谢！

wuhaimiao

万分感谢！