|
|
真的很谢谢你！
如开篇程序所言在MATLAB10版本中对行列式总体求导所得的解析式也是能够较快求出，但在赋值时迟迟没有进展；但分开求导，哪怕所得解析式仍汇合在一起 ，最后对总解析式赋值，却仍能较快求出数值。不知何故？？？
两种方法对最后赋值的影响难道机理上有什么玄机？
clear all
format long
syms u v
m=1;
s=1+0.1;
a=0.98;
b=0.99;
n1=1.48;
n2=-a*n1;
n3=b*n1;
r=(n1^2-n2^2)/(n1^2-n3^2);
w1=v^2*r-u^2;
w=v^2-u^2;
ne=sqrt((1-u^2/v^2)*(n1^2-n3^2)+n3^2);
F1=-besselj(m+1,u)+m./u.*besselj(m,u);
F=F1./(u.*besselj(m,u));
M1=besseli(m+1,w1)+m./w1.*besseli(m,w1);
M=M1./(w1.*besseli(m,w1));
N1=-besselk(m+1,w1)+m./w1.*besselk(m,w1);
N=N1./(w1.*besselk(m,w1));
Q=m*(1./u.^2+1./w1.^2);
R=m*(1./w1.^2-1./w.^2)/s;
S1=besseli(m+1,w1*s)+m./w1./s.*besseli(m,w1*s);
S=S1./(w1.*besseli(m,w1*s));
T1=-besselk(m+1,w*s)+m./w./s.*besselk(m,w*s);
T=T1./(w.*besselk(m,w*s));
X1=-besselk(m+1,w1*s)+m./w1./s.*besselk(m,w1*s);
X=X1./(w1.*besselk(m,w1*s));
I1m=besseli(m,w1);
K1m=besselk(m,w1);
I2m=besseli(m,w1*s);
K2m=besselk(m,w1*s);
x1=Q.*ne.*I1m;
x2=Q.*ne.*K1m;
x3=I1m.*(F-M);
x4=K1m.*(F-N);
x5=R*ne.*I2m;
x6=R.*ne.*K2m;
x7=I2m.*(-S-T);
x8=K2m.*(-X-T);
x9=I1m.*(n1^2.*F-n2^2.*M);
x10=K1m.*(n1^2.*F-n2^2.*N);
x11=Q.*ne.*I1m;
x12=Q.*ne.*K1m;
x13=I2m.*(-S.*n2^2-n3^2.*T);
x14=K2m.*(-n2^2.*X-n3^2.*T);
x15=R.*ne.*I2m;
x16=R.*ne.*K2m;
% f=det([x1,x2,x3,x4;
% x5,x6,x7,x8;
% x9,x10,x11,x12;
% x13,x14,x15,x16]);
%%PSL@CSU
%%QQ:547423688
%%Email:anyuezhiji@qq.com
%%Edit @ 2010.4.17
%拆开算
detf1=x1*x6*x11*x16;
detf2=-x1*x6*x12*x15;
detf3=-x1*x10*x7*x16;
detf4=x1*x10*x8*x15;
detf5=x1*x14*x7*x12;
detf6=-x1*x14*x8*x11;
detf7=-x5*x2*x11*x16;
detf8=x5*x2*x12*x15;
detf9=x5*x10*x3*x16;
detf10=-x5*x10*x4*x15;
detf11=-x5*x14*x3*x12;
detf12=x5*x14*x4*x11;
detf13=x9*x2*x7*x16;
detf14=-x9*x2*x8*x15;
detf15=-x9*x6*x3*x16;
detf16=x9*x6*x4*x15;
detf17=x9*x14*x3*x8;
detf18=-x9*x14*x4*x7;
detf19=-x13*x2*x7*x12;
detf20=x13*x2*x8*x11;
detf21=x13*x6*x3*x12;
detf22=-x13*x6*x4*x11;
detf23=-x13*x10*x3*x8;
detf24=x13*x10*x4*x7;
for i=1:24
eval(['fv',num2str(i),'=diff(detf',num2str(i),',v);']);
end
for i=1:24
eval(['fu',num2str(i),'=diff(detf',num2str(i),',u);']);
end
%这里开始读入数据
% v=6.4414247e-001
%u=6.4091366e-001
%q=load('e:\0.1.txt');
%for i=1:length(q(:,1))
% u=q(i,3)
% v=q(i,1)
%下面开始求值
ffv=0;
for i=1:24
eval(['ffv=ffv+eval(fv',num2str(i),');'])
end
ffu=0;
for i=1:24
eval(['ffu=ffu+eval(fu',num2str(i),');'])
end
ff1=-ffv/ffu;
ff2=1+(u/v)^2*(1-2*v/u*ff1);
q=load('e:\0.1.txt');
for ii=1:length(q(:,1))
v=q(ii,1)
u=q(ii,3)
ff=eval(ff2)
hold on
plot(v,ff)
end |
|
回复此楼