声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1078|回复: 5

[图像处理] 图形为什么没有原文献中的光滑

[复制链接]
发表于 2008-3-12 23:52 | 显示全部楼层 |阅读模式

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

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

x
这是我的程序,大家帮忙看看,得到图形为什么没有原文献中的光滑,是数据精度的问题,还是程序写得太笨了,谢谢!!

clear all;clc;close all;
global Tc;
sz1 = [];
sz2 = [];
sz3 = [];
sz4 = [];
sz5 = [];
Tc_p = 0.01:0.01:2.5;
n0= [0.2 0.2 0.2 0.2 0.2];
    Tc = Tc_p(1);
    options=optimset('display','on');
    n(1,:)=fsolve(@myfun,n0,options);
    for i = 2:length(Tc_p)
        Tc = Tc_p(i);
        n(i,:)=fsolve(@myfun,n(i-1,:),options);
    end
    sz1 = [sz1 n(:,1)];
    sz2 = [sz2 n(:,2)];
    sz3 = [sz3 n(:,3)];
    sz4 = [sz4 n(:,4)];
    sz5 = [sz5 n(:,5)];

syms k1 k2 k3 k4 k5
syms x1 x2 x3 x4 x5
syms T
syms a b c d e
X=[];
js=1.0;
ja=1.0;
je=1.0;
qs=1.0;
q=3.5;
k=1.0;
f1=x1/(2*sqrt(qs^2+x1^2))*tanh(sqrt(qs^2+x1^2)/(2*k*T));
f2=x2/(2*sqrt(q^2+x2^2))*tanh(sqrt(q^2+x2^2)/(2*k*T));
f3=x3/(2*sqrt(q^2+x3^2))*tanh(sqrt(q^2+x3^2)/(2*k*T));
f4=x4/(2*sqrt(q^2+x4^2))*tanh(sqrt(q^2+x4^2)/(2*k*T));
f5=x5/(2*sqrt(q^2+x5^2))*tanh(sqrt(q^2+x5^2)/(2*k*T));
d1=diff(f1,x1);
d2=diff(f2,x2);
d3=diff(f3,x3);
d4=diff(f4,x4);
d5=diff(f5,x5);

for i = 1:length(Tc_p)
   
    R1=sz1(i);R2=sz2(i);R3=sz3(i);R4=sz4(i);R5=sz5(i);
   d11=subs(d1,T,Tc_p(i));
   d22=subs(d2,T,Tc_p(i));
   d33=subs(d3,T,Tc_p(i));
   d44=subs(d4,T,Tc_p(i));
   d55=subs(d5,T,Tc_p(i));
   
   d111=subs(d11,x1,4*js*R1+je*R2);
   d222=subs(d22,x2,4*ja*R2+je*R1+je*R3);
   d333=subs(d33,x3,4*ja*R3+je*R2+je*R4);
   d444=subs(d44,x4,4*ja*R4+je*R3+je*R5);
   d555=subs(d55,x5,4*ja*R5+je*R4+je*R5);

   y1=(1/d111-4*js)*k1-je*k2-1;
   y2=-je*k1+(1/d222-4*ja)*k2-je*k3-1;
   y3=-je*k2+(1/d333-4*ja)*k3-je*k4-1;
   y4=-je*k3+(1/d444-4*ja)*k4-je*k5-1;
   y5=-je*k4+(1/d555-4*ja)*k5-je*k5-1;

   [a,b,c,d,e]=solve(y1,y2,y3,y4,y5,'k1','k2','k3','k4','k5');
   a=double(a);
   b=double(b);
   c=double(c);
   d=double(d);
   e=double(e);

   X(i)=10/(2/(1+a)+2/(1+b)+2/(1+c)+2/(1+d)+2/(1+e))-1;

end
plot(Tc_p,X);


%%%%%%%%%%%%下面是调用函数myfun.m%%%%%%%%%%%
function y=myfun(n)
global Tc
syms sz1 sz2 sz3 sz4 sz5
js=1.0;
ja=1.0;
je=1.0;
qs=1.0;
q=3.5;
k=1.0;
sz1=n(1);
sz2=n(2);
sz3=n(3);
sz4=n(4);
sz5=n(5);
t1=sqrt(qs^2+(4*js*sz1+je*sz2)^2);
t2=sqrt(q^2+(4*ja*sz2+je*sz1+je*sz3)^2);
t3=sqrt(q^2+(4*ja*sz3+je*sz2+je*sz4)^2);
t4=sqrt(q^2+(4*ja*sz4+je*sz3+je*sz5)^2);
t5=sqrt(q^2+(4*ja*sz5+je*sz4+je*sz5)^2);
y(1)=(4*js*sz1+je*sz2)/(2*t1)*tanh(t1/(2*k*Tc))-sz1;
y(2)=(4*ja*sz2+je*sz1+je*sz3)/(2*t2)*tanh(t2/(2*k*Tc))-sz2;
y(3)=(4*ja*sz3+je*sz2+je*sz4)/(2*t3)*tanh(t3/(2*k*Tc))-sz3;
y(4)=(4*ja*sz4+je*sz3+je*sz5)/(2*t4)*tanh(t4/(2*k*Tc))-sz4;
y(5)=(4*ja*sz5+je*sz4+je*sz5)/(2*t3)*tanh(t5/(2*k*Tc))-sz5;


运行程序后,得到图形在峰的位置和高度都是和原文献一样的,但是在横坐标 x=1 附近为什么不很光滑,而原文献是很光滑的,是最后a,b,c,d,e的数据精度的问题吗?还是程序写的太笨了?

麻烦大家了,谢谢啊谢谢啊!!
回复
分享到:

使用道具 举报

发表于 2008-3-13 09:25 | 显示全部楼层
程序那么长,恐怕只有你自己知道了
发表于 2008-3-13 09:35 | 显示全部楼层
这些问题建议还是自己解决吧,你觉得哪里有问题,就显示哪个范围的数据,或者设置断点调试一下
 楼主| 发表于 2008-3-13 10:09 | 显示全部楼层

回复 3楼 的帖子

谢谢
 楼主| 发表于 2008-3-24 15:39 | 显示全部楼层

再次求助:图形为什么没有原文献的光滑?

下面是我的程序,程序是用fsolve解了一个非线性方程组,然后作图,可是得到的图形没有原文献的光滑?想了很长时间,不知道为什么?两个图形在附件中。
clear all;clc;close all;
global Tc;
global je
sz1 = [];
sz2 = [];
je_p = [0.05 0.5 1.0];
Tc_p = 0.0:0.01:2.0;
n0= [0.3 0.3];
for je = je_p
    Tc = Tc_p(1);
    options=optimset('display','on');
    n(1,:)=fsolve(@myfun7,n0,options);
    for i = 2:length(Tc_p)
        Tc = Tc_p(i);
        n(i,:)=fsolve(@myfun7,n(i-1,:),options);
    end
    sz1 = [sz1 n(:,1)];
    sz2 = [sz2 n(:,2)];
end
plot(Tc_p,sz1);
hold on
plot(Tc_p,sz2);
##################下面是myfun7.m################
function y=myfun7(n)
global Tc
global je
syms sz1 sz2
js=1.3;
ja=1.0;
qs=2.0;
q=1.0;
k=1.0;
sz1=n(1);
sz2=n(2);
Hs=sqrt(qs^2+(4*js*sz1+je*sz2)^2);
H=sqrt(q^2+(4*ja*sz2+2*je*sz1)^2);
y(1)=(4*js*sz1+je*sz2)/(2*Hs)*tanh(Hs/(2*k*Tc))-sz1;
y(2)=(4*ja*sz2+2*je*sz1)/(2*H)*tanh(H/(2*k*Tc))-sz2;

Doc1.doc

42.5 KB, 下载次数: 10

 楼主| 发表于 2008-3-24 15:41 | 显示全部楼层

回复 楼主 的帖子

考虑好长时间了,感觉离“事实的真相'不远了,就是和原图有些差别,请有这方面经验的帮小弟看看,谢谢!!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 03:23 , Processed in 0.067878 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表