请教！请教！精度问题

stu000 · 发表于 2006-4-13 14:00

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

用2分法计算时，同一个程序什么都不改动，居然2次搜出来的根不一样，怀疑是精度问题，哪位好心的朋友能给解答一下呀，该怎么改呢！！！？？？小女子万分感谢！！！[em06]

aspen · 发表于 2006-4-13 17:54

试试vpa

ilfl · 发表于 2006-4-14 15:51

使用format long可将精度提高

happy · 发表于 2006-4-14 15:53

<DIV class=quote>以下是引用ilfl在2006-4-14 15:51:41的发言： 使用format long可将精度提高 </DIV>
 这个只是显示精度

xjzhang · 发表于 2006-4-14 20:35

能不能换一种求根方法：solve ,roots ,fzero

happy · 发表于 2006-4-15 08:39

<DIV class=quote>以下是引用xjzhang在2006-4-14 20:35:00的发言： 
能不能换一种求根方法：solve ,roots ,fzero</DIV>
 什么意思？

stu000 · 发表于 2006-4-16 14:29

谢谢各位好心的哥哥！！！不好意思，回复晚了呀：）小女子把程序贴出来，哪位好心的哥哥帮我看看呀！ global x ; xhr=zeros(200); vel=zeros(200,100); xtol=1.0e-6; xhs=0.25; xhe=0.25; nx=1; cs=0.90 ce=0.905 nc=20; [xhr,vel]=velsearch(nx,xhs,xhe,nc,cs,ce,xtol); % Plot the velocity verses the length ratio. plot(xhr,vel,'k.') axis([xhs xhe cs ce]) a=vel

stu000 · 发表于 2006-4-16 14:33

function fp=eqn_w3(v,xh); % xh=0.25; % v=0.90;
xhh=0.25; roub=19300; rou=3879; c11b=18.6*10^10; c44b=4.2*10^10; c11=6.13*10^10; c44=2.18*10^10; lambda=c11-2*c44; lambdab=c11b-2*c44b; mu=c44; mub=c44b;
l2m = lambda + 2 * mu;
l2mb = lambdab + 2 * mub; % 1.不带电极层区域 hhhh=10; H=0.002;
 beta1=0.83797530; beta2=0.40325621;
beta(1)=beta1; beta(2)=beta2; Beta1=0.84081852; Beta2=0.41956813; Beta(1)=Beta1; Beta(2)=Beta2;

% Integration of the thickness expansion functions for large depth h fact1 = 1 - 1 / exp(2 * pi * (beta1 + beta1) * hhhh); fact2 = 1 - 1 / exp(2 * pi * (beta2 + beta2) * hhhh); fact12= 1 - 1 / exp(2 * pi * (beta1 + beta2) * hhhh); f11 = fact1 / (beta1 + beta1); f12 = fact12 / (beta1 + beta2); f21 = f12; f22 = fact2 / (beta2 + beta2);
% Define the frequency matrix A(4, 4) syms k; A(1, 1) =f11 * (l2m * k^2 + mu * beta1^2 - mu * v^2 ); A(1, 2) =f11 * beta1 * (lambda - mu) * k; A(1, 3) =f21 * (l2m * k^2 + mu * beta1 * beta2 - mu * v^2); A(1, 4) =f21 * (lambda * beta2 - mu * beta1) * k; A(2, 1) =f11 * beta1 * (lambda - mu) * k; A(2, 2) =f11 * (mu * k^2 + l2m * beta1^2 - mu * v^2 ); A(2, 3) =f21 * (beta1 * lambda - beta2 * mu) * k; A(2, 4) =f21 * (mu * k^2 + l2m * beta1 * beta2 - mu * v^2); A(3, 1) =f21 * (l2m * k^2 + mu * beta1 * beta2 - mu * v^2); A(3, 2) =f21 * (beta1 * lambda - beta2 * mu) * k; A(3, 3) =f22 * (l2m * k^2 + mu * beta2^2 - mu * v^2); A(3, 4) =f22 * beta2 * (lambda - mu) * k; A(4, 1) =f21 * (lambda * beta2 - mu * beta1) * k; A(4, 2) =f21 * (mu * k^2 + l2m * beta1 * beta2 - mu * v^2); A(4, 3) =f22 * beta2 * (lambda - mu) * k; A(4, 4) =f22 * (mu * k^2 + l2m * beta2^2 - mu * v^2);
y=det(A);
 kk=solve(y,k);
 G(1)=numeric(real(kk(1))); G(2)=numeric(real(kk(3))); G(3)=numeric(imag(kk(5))*i); G(4)=numeric(imag(kk(7))*i); 
 
for n=1:4, d(1,1)=f11 * (l2m * G(n)^2 + mu * beta1^2 - mu * v^2 ); d(1,2)=f11 * beta1 * (lambda - mu) * G(n) ; d(1,3)=f21 * (l2m * G(n)^2 + mu * beta1 * beta2 - mu * v^2); d(2,1)=f11 * beta1 * (lambda - mu) * G(n) ; d(2,2)=f11 * (mu * G(n)^2 + l2m * beta1^2 - mu * v^2 ); d(2,3)=f21 * (beta1 * lambda - beta2 * mu) * G(n); d(3,1)=f21 * (l2m * G(n)^2 + mu * beta1 * beta2 - mu * v^2); d(3,2)=f21 * (beta1 * lambda - beta2 * mu) * G(n); d(3,3)=f22 * (l2m * G(n)^2 + mu * beta2^2 - mu * v^2); h(1,1)=-f21 * (lambda * beta2 - mu * beta1) * G(n);
 h(2,1)=-f21 * (mu * G(n)^2 + l2m * beta1 * beta2 - mu * v^2);
 h(3,1)=-f22 * beta2 * (lambda - mu) * G(n);
al(:,n)=d\h;
end;
 % new lame constants for m=1:2 for n = 1 : 2 lambdabb(m,n) = lambda * (1+2*H*(Beta(m)+Beta(n))*lambdab/lambda); mubb(m,n)=mu*(1+2*H*(Beta(m)+Beta(n))*mub/mu); l2mbb(m,n)=l2m*(1+2*H*(Beta(m)+Beta(n))*l2mb/l2m); rour(m,n)=(1+2*(Beta(m)+Beta(n))*H*roub/rou); end; end;
 % 基体厚度hhhh 电极厚度H （归一化后） % Integration of the thickness expansion functions for large depth h Fact1 = 1 - 1 / exp(2 * pi * (Beta1 + Beta1) * hhhh); Fact2 = 1 - 1 / exp(2 * pi * (Beta2 + Beta2) * hhhh); Fact12= 1 - 1 / exp(2 * pi * (Beta1 + Beta2) * hhhh); F11 = Fact1 / (Beta1 + Beta1); F12 = Fact12 / (Beta1 + Beta2); F21 = F12; F22 = Fact2 / (Beta2 + Beta2);
% Define the frequency matrix A(4, 4) syms K; AA(1, 1) = F11 * (l2mbb(1,1) * K^2 + mubb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); AA(1, 2) = + F11 * Beta1 * (lambdabb(1,1) - mubb(1,1)) * K ; AA(1, 3) = F21 * (l2mbb(1,2) * K^2 + mubb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); AA(1, 4) = + F21 * (lambdabb(1,2) * Beta2 - mubb(1,2) * Beta1) * K; AA(2, 2) = F11 * (mubb(1,1) * K^2 + l2mbb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); AA(2, 3) = + F21 * (Beta1 * lambdabb(1,2) - Beta2 * mubb(1,2)) * K; AA(2, 4) = F21 * (mubb(1,2) * K^2 + l2mbb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); AA(3, 3) = F22 * (l2mbb(2,2) * K^2 + mubb(2,2) * Beta2^2 - mu * v^2 * rour(2,2)); AA(3, 4) = + F22 * Beta2 * (lambdabb(2,2) - mubb(2,2)) * K; AA(4, 4) = F22 * (mubb(2,2) * K^2 + l2mbb(2,2) * Beta2^2 - mu * v^2 * rour(2,2)); AA(2, 1) = AA(1, 2); AA(3, 1) = AA(1, 3); AA(4, 1) = AA(1, 4); AA(4, 3) = AA(3, 4); AA(3, 2) = AA(2, 3); AA(4, 2) = AA(2, 4); yy = det(AA); % 并求振幅比AL(m,n) KK=solve(yy,K); GG(1)=numeric(real(KK(1))); GG(2)=numeric(real(KK(3))); GG(3)=numeric(imag(KK(5))*i); GG(4)=numeric(imag(KK(7))*i);
for n=1:4, dd(1,1)=F11 * (l2mbb(1,1) * GG(n)^2 + mubb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); dd(1,2)=F11 * Beta1 * (lambdabb(1,1) - mubb(1,1)) * GG(n) ; dd(1,3)=F21 * (l2mbb(1,2) * GG(n)^2 + mubb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); dd(2,1)=F11 * Beta1 * (lambdabb(1,1) - mubb(1,1)) * GG(n); dd(2,2)=F11 * (mubb(1,1) * GG(n)^2 + l2mbb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); dd(2,3)=F21 * (Beta1 * lambdabb(1,2) - Beta2 * mubb(1,2)) * GG(n); dd(3,1)=F21 * (l2mbb(1,2) * GG(n)^2 + mubb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); dd(3,2)=F21 * (Beta1 * lambdabb(1,2) - Beta2 * mubb(1,2)) * GG(n); dd(3,3)=F22 * (l2mbb(2,2) * GG(n)^2 + mubb(2,2) * Beta2^2 - mu * v^2 * rour(2,2)); hh(1,1)=-F21 * (lambdabb(1,2) * Beta2 - mubb(1,2) * Beta1) * GG(n);
 hh(2,1)=-F21 * (mubb(1,2) * GG(n)^2 + l2mbb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2));
 hh(3,1)=-F22 * Beta2 * (lambdabb(2,2) - mubb(2,2)) * GG(n);
 Al(:,n)=dd\hh;
 end;

 for n=1:4 ddd(1,1)=F11 * (l2mbb(1,1) * GG(n)^2 + mubb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); ddd(1,2)=-F11 * Beta1 * (lambdabb(1,1) - mubb(1,1)) * GG(n) ; ddd(1,3)=F21 * (l2mbb(1,2) * GG(n)^2 + mubb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); ddd(2,1)=-F11 * Beta1 * (lambdabb(1,1) - mubb(1,1)) * GG(n); ddd(2,2)=F11 * (mubb(1,1) * GG(n)^2 + l2mbb(1,1) * Beta1^2 - mu * v^2 * rour(1,1)); ddd(2,3)=-F21 * (Beta1 * lambdabb(1,2) - Beta2 * mubb(1,2)) * GG(n); ddd(3,1)=F21 * (l2mbb(1,2) * GG(n)^2 + mubb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2)); ddd(3,2)=-F21 * (Beta1 * lambdabb(1,2) - Beta2 * mubb(1,2)) * GG(n); ddd(3,3)=F22 * (l2mbb(2,2) * GG(n)^2 + mubb(2,2) * Beta2^2 - mu * v^2 * rour(2,2)); hhh(1,1)=F21 * (lambdabb(1,2) * Beta2 - mubb(1,2) * Beta1) * GG(n);
 hhh(2,1)=-F21 * (mubb(1,2) * GG(n)^2 + l2mbb(1,2) * Beta1 * Beta2 - mu * v^2 * rour(1,2));
 hhh(3,1)=F22 * Beta2 * (lambdabb(2,2) - mubb(2,2)) * GG(n);
 Bl(:,n)=ddd\hhh; end;
 % 3.验证边界条件 for m=1:4 zterm=pi*GG(m)*xhh; ztermm=pi*G(m)*xhh; ztermmm=pi*G(m)*(xh+xhh); [sinQQ, cosQQ] = trigimag(zterm); [sinU,cosU]=trigimag(ztermm); [sinQ,cosQ]=trigimag(ztermmm);
T(1,m)=((F11*l2mbb(1,1)*Al(1,m)+F21*l2mbb(1,2)*Al(3,m))*GG(m)+(F11*lambdabb(1,1)*Beta1*Al(2,m)+F21*lambdabb(1,2)*Beta2))*cosQQ; T(1,m+4)=(-(F11*l2mbb(1,1)*Bl(1,m)+F21*l2mbb(1,2)*Bl(3,m))*GG(m)+(F11*lambdabb(1,1)*Beta1*Bl(2,m)+F21*lambdabb(1,2)*Beta2))*sinQQ; T(1,m+8)=-(l2m * (f11 * al(1, m) + f12 * al(3, m))*G(m) + lambda * (f11 * beta1 * al(2, m) + f12 * beta2)) * cosU;

T(2,m)=((-F11*mubb(1,1)*Al(2,m)-F21*mubb(1,2))*GG(m)+F11*mubb(1,1)*Beta1*Al(1,m)+F21*mubb(1,2)*Beta2*Al(3,m))*sinQQ; T(2,m+4)=((F11*mubb(1,1)*Bl(2,m)+F21*mubb(1,2))*GG(m)+F11*mubb(1,1)*Beta1*Bl(1,m)+F21*mubb(1,2)*Beta2*Bl(3,m))*cosQQ; T(2,m+8)=-mu*((-f11 * al(2, m)-f21)*G(m) + f11 * beta1 * al(1, m) + f21 * beta2* al(3, m)) * sinU; 
 T(3,m)=((F21*l2mbb(1,2)*Al(1,m)+F22*l2mbb(2,2)*Al(3,m))*GG(m)+(F21*lambdabb(1,2)*Beta1*Al(2,m)+F22*lambdabb(2,2)*Beta2))*cosQQ; T(3,m+4)=(-(F21*l2mbb(1,2)*Bl(1,m)+F22*l2mbb(2,2)*Bl(3,m))*GG(m)+(F21*lambdabb(1,2)*Beta1*Bl(2,m)+F22*lambdabb(2,2)*Beta2))*sinQQ; T(3,m+8)=-(l2m *(f21 * al(1,m)+f22*al(3,m))*G(m)+lambda*(f21 * beta1 * al(2, m) + f22 * beta2)) * cosU;
 
T(4,m)=((-F21*mubb(1,2)*Al(2,m)-F22*mubb(2,2))*GG(m)+F21*mubb(1,2)*Beta1*Al(1,m)+F22*mubb(2,2)*Beta2*Al(3,m))*sinQQ; T(4,m+4)=((F21*mubb(1,2)*Bl(2,m)+F22*mubb(2,2))*GG(m)+F21*mubb(1,2)*Beta1*Bl(1,m)+F22*mubb(2,2)*Beta2*Bl(3,m))*cosQQ; T(4,m+8)=-mu*((-f21 * al(2, m)-f22)*G(m) + f21 * beta1 * al(1, m) + f22 * beta2* al(3, m)) * sinU;
 T(5,m)=Al(1,m)*sinQQ; T(5,m+4)=Bl(1,m)*cosQQ; T(5,m+8)=-al(1,m)*sinU;
 T(6,m)=Al(2,m)*cosQQ; T(6,m+4)=Bl(2,m)*sinQQ; T(6,m+8)=-al(2,m)*cosU;

T(7,m)=Al(3,m)*sinQQ; T(7,m+4)=Bl(3,m)*cosQQ; T(7,m+8)=-al(3,m)*sinU;
 T(8,m)=cosQQ; T(8,m+4)=sinQQ; T(8,m+8)=-cosU;
T(9,m)=0; T(9,m+4)=0; T(9,m+8)=-al(1,m)*sinQ;
T(10,m)=0; T(10,m+4)=0; T(10,m+8)=-al(3,m)*sinQ;
T(11,m)=0; T(11,m+4)=0; T(11,m+8)=-mu*((-f11 * al(2, m)-f21)*G(m) + f11 * beta1 * al(1, m) + f21 * beta2* al(3, m)) * sinQ;
 T(12,m)=0; T(12,m+4)=0; T(12,m+8)=-mu*((-f21 * al(2, m)-f22)*G(m) + f21 * beta1 * al(1, m) + f22 * beta2* al(3, m)) * sinQ; end 
fp=numeric(det(T))

stu000 · 发表于 2006-4-16 14:34

function [xhr,vel]=velsearch(nx,xhs,xhe,nc,cs,ce,xtol) %global x c c6; % set up the biscection search procedure. % xh0=20.0; %xtol=1.0e-10 ir=0; % veldat=fopen('vel.dat','w'); % nx=1; for i=1:nx, v=xhs+(i-1)*(xhe-xhs)/nx; % consider the at-cut of quartz (32.15 Y-X) % phi=pi*xh/180; % [c,c6]=elastranyx(phi); %nc=50 %v=40 % cs=0.9;ce=1.2; x=cs ds=(ce-cs)/nc; dso=ds; fo=(eqn_w3(x,v)); if abs(real(fo))<abs(imag(fo)); fo=imag(fo); else fo=real(fo); end; %[fo,beta]=wwavenums1(xh); ir=0; for j=1:1000 x=x+ds; if x>ce break;end; fp=(eqn_w3(x,v)); if abs(real(fp))<abs(imag(fp)); fp=imag(fp); else fp=real(fp); end; %[fp,beta]=wwavenums1(xh); %fprintf(JJJ=%4i XXX=%13.10f OOO=%10.4e FFF=%10.4e DDD=%10.4e\n',... % j,x,fo,fp,ds) if (fo*fp>0) fo=fp; %elseif (ds<xtol &abs(fp)<xtol) elseif (ds<xtol) % 'XH= ',xh,'Root= ',x %fprintf('IRRR= %4i XH= %12.8f Root=%15.12f \n',ir,xh,x); %fprintf('Beta(1)=%15.12f\n real(2)=%15.12f\n Imag(2)=%15.12f\n',... % beta(1),real(beta(2)),abs(imag(beta(2)))); %fp,fo %fprintf(veldat,' %10.4f %13.10f\n',xh,x); %break; ir=ir+1; xhr(i)=v vel(i,ir)=x save velxhr.mat xhr vel; fo=fp; ds=dso; %break; else %fprint('JJJ= %4i XXX=%13.10f OOO=%10.4e FFF=%10.4e DDD=%10.4e\n',... %j,x,fo,fp,ds) x=x-ds; ds=ds/2; end; end; end; %fclose(veldat); end end

stu000 · 发表于 2006-4-16 14:34

function [sinl, cosl] = trigimag(zterm)
cosl = 0; sinl = 0;
% No scaling is needed if the argument is real. if (imag(zterm) == 0) sinl = sin(zterm); cosl = cos(zterm); % Scale the sine and cosine functions in case of imagnary argument elseif (imag(zterm) > 0) %exp(imag(zterm)) chl = 0.5 * (1 + exp( - 2 * imag(zterm))); shl = 0.5 * (1 - exp( - 2 * imag(zterm))); cosl = cos(real(zterm)) * chl - i * sin(real(zterm)) * shl; sinl = sin(real(zterm)) * chl + i * cos(real(zterm)) * shl; else %exp( - imag(zterm)) chl = 0.5 * ( 1 + exp(2 * imag(zterm))); shl = 0.5 * (- 1 + exp(2 * imag(zterm))); cosl = cos(real(zterm)) * chl - i * sin(real(zterm)) * shl; sinl = sin(real(zterm)) * chl + i * cos(real(zterm)) * shl; end; clear chl shl

stu000 · 发表于 2006-4-16 14:40

vpa命令也都用了，还是不行 现在就是仅仅在xh=0.25这个点在0.90---0.905之间来搜根，一般会搜出4，5个根左右，但是我同一个程序还是同一个点，搜索范围也一样，什么都不变，，前后两次运行同一个程序，居然搜索出来的根跟前一次不一样，超级郁闷，应该是计算程序的精度问题，但是小女子实在找不出错误在哪里，麻烦哪位好心的大哥一定费心帮我看看，我现在搞的是毕业论文，咳，要是弄不出来就毕业不了啊，郁闷：（麻烦各位了，定有报答：）

[此贴子已经被作者于2006-4-16 15:32:06编辑过]

stu000 · 发表于 2006-4-17 10:19

顶。。。。。。。。最后算出的det(T)值很大啊，有10的60次方，我想把程序中比较大的数都除以一个十的几次方什么的。 按理说我前面把精度提高到了150，应该可以算出相同的结果啊，结果还是不行 郁闷。。。。。。。。。。。

happy · 发表于 2006-4-17 11:09

个人认为，你这个问题不是精度的问题，还是算法或者程序的问题 即使精度不够高，但是两次计算你没有改变任何东西且程序中不存在随机性因素在里边 那迭代也不会出现两个不同的解。
程序比较长，没时间帮你检查，我想这个工作应该你自己来完成
所以建议还是好好查查程序，看程序中是否存在随机性因素，可以用一个简单的模型验证一下程序的正确性
另外根据你的说法，这个方程应该是个病态方程，可以考虑除以个什么数，当然也可以换个算法，比如超松弛迭代之类的。

stu000 · 发表于 2006-4-19 13:46

谢谢happy老师，我再仔细检查一下！！！

Hamlet · 发表于 2006-4-19 16:59

有点晕呐

账号		自动登录	找回密码
密码			我要加入

请教！请教！精度问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复：（stu000）请教！请教！精度问题

回复：（ilfl）使用format long可将精度提高

回复：（xjzhang）能不能换一种求根方法：solve ,root...

回复：（stu000）请教！请教！精度问题

...

浏览过的版块