| 查看: 2153 | 回复: 15 | ||
[求助]
拟合三维曲线 已有1人参与
|
|
最近学习拟合三维曲线，之前是小白。毕业师哥给的程序，数据是从ＡＮＳＵＳ温度场提出的温度为９２０°Ｃ的坐标，现将其拟合处三维曲面。图片和文件上传老是失败 有大神留言加个ＱＱ我发你看看吧，万分感谢！ 我自己拟合了几次 程序一直出错，求大神帮看看吧，最终完美的图片我附上，是不是方法错了？ >> data=load('920.txt'); y=data(:,1); x=data(:,3); z=data(:,2); >> [X,Y]=meshgrid(-0.25:0.001:0.1,-0.015:0.0005:0.015) Z=griddata(x,y,z,X,Y) mesh(X,Y,Z) plot3(x,y,z,'*') pcolor(X,Y,Z);shading interp%伪彩色图 figure,contourf(X,Y,Z) %等高线图 figure,surf(X,Y,Z);%三维曲面 data=load('920.txt'); x=data(:,1:2); y=data(:,3); beat0=[1,1,1,1]'; [beat,r,j]=nlinfit(x,y,@myfunn,beat0) data=load('9201.txt'); b=data(:,1); a=data(:,3); c=data(:,2); -c*abs(sqrt(1-(x-d)^2/a^2-y^2/b^2)) data=load('TXYZ.txt'); b=data(:,1); a=data(:,3); c=data(:,2); d=data(:,4); f=@(q,x)1.87.*q(5).*exp(-3.*((x(:,1)-q(4)).^2./q(1).^2+x(:,2).^2./q(2).^2+x(:,3).^2./q(3).^2))./(q(1).*q(2).*q(3))-x(:,4) f=@(q,x)exp(abs(sqrt((x(:,1)-q(1)).^2+x(:,2).^2+x(:,3).^2))-(x(:,1)-q(1)))-q(2).*abs(sqrt((x(:,1)-q(1)).^2+x(:,2).^2+x(:,3).^2)) f=@(q,x)exp(abs(sqrt(x(:,1).^2+x(:,2).^2+x(:,3).^2))-x(:,1))-q(1).*abs(sqrt(x(:,1).^2+x(:,2).^2+x(:,3).^2)) >> x=[a b c]; >> q0=[1 1 1 1 1]; q0=[0.2 0.00637 0.0038 -0.12 0.001]; q=nlinfit(x,zeros(size(x,1),1),f,q0); disp(num2str(q)); H A B分别为 3.388 0.39101 -1.5654 3.388 f=@(x,y,z)1.87.*0.11689.*exp(-3.*((x+0.054754).^2./0.53846.^2+y.^2./0.042682.^2+z.^2./0.014206.^2))./(0.53846.*0.042682.*0.014206)-550 f=@(x,y,z)exp(abs(sqrt(x.^2+y.^2+z.^2))-x)-20.7989.*abs(sqrt(x.^2+y.^2+z.^2)) f=@(x,y,z)exp(-1.5654.*(x-sqrt(x.^2+y.^2+z.^2))./(2.*0.39101))-3.388.*3.388.*sqrt(x.^2+y.^2+z.^2) >> [x,y,z] = meshgrid(-0.4:0.005:0.4,-0.02:0.0004:0.02,-0.005:0.0001:0); >> v = f(x,y,z); h = patch(isosurface(x,y,z,v,0)); >> isonormals(x,y,z,v,h) set(h,'FaceColor','r','EdgeColor','none'); xlabel('x');ylabel('y');zlabel('z'); alpha(1) grid on; view([1,1,1]); axis equal; camlight; lighting gouraud >> >> hold on -c*abs(sqrt(1-(x-d)^2/a^2-y^2/b^2)) f=@(x,y,z)z+0.0032*abs(sqrt(1-x^2/0.08917^2-y^2/0.0032^2)) [x,y] = meshgrid(0:0.002:0.08917,-0.00527:0.0002:0.00527); f=@(x,y,z)(x+0.0295).^2/0.08917.^2+y.^2/0.00527.^2+z.^2/0.0032.^2-1 [x,y,z] = meshgrid(-0.04:0.002:0.062,-0.005:0.0002:0.005,-0.0032:0.0001:0); v = f(x,y,z); h = patch(isosurface(x,y,z,v,0)); >> isonormals(x,y,z,v,h) set(h,'FaceColor','r','EdgeColor','none'); xlabel('x');ylabel('y');zlabel('z'); alpha(1) grid on; view([1,1,1]); axis equal; camlight; lighting gouraud x=linspace(-0.04,0.062,50); y=linspace(-0.005,0.005,50); [xx,yy]=meshgrid(x,y); zz=-0.0032.*abs(sqrt(1-xx.^2/0.08917.^2-yy.^2/0.0032.^2)) i=find(yy<=-3&yy>=3); z1=zz;z1(i)=nan; surf(xx,yy,z1) x=linspace(-0.04,0.062,50); y=linspace(-0.005,0.005,50); [xx,yy]=meshgrid(x,y); zz=-0.0032.*abs(sqrt(1-xx.^2/0.08917.^2-yy.^2/0.0032.^2)) surf(xx,yy,zz);shading interp; pause i=find(yy<=-0.003|yy>=0.003); aa=zz;aa(i)=NaN; g=find(yy>=0.003); aa=zz;aa(g)=NaN; surf(xx,yy,aa);shading interp i=find(yy.^2>=0.000009); aa=zz;aa(i)=NaN; surf(xx,yy,aa);shading interp  Y-Z方向.jpg |
回复此楼» 猜你喜欢
93BebMhtakh前后11位开头都是大写
已经有5人回复
-
欢迎发来filecode的Mz6后的代码验证其规律
已经有86人回复
-
感觉是下周放榜了
已经有11人回复
-
filecode=后面第一个是大写字母
已经有10人回复
-
时间戳变了，能看出什么问题？
已经有4人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有15人回复
-
今天系统多次维护，明天很可能放榜！
已经有9人回复
-
时间戳又变了8-15
已经有25人回复
-
2027广东省杰青
已经有6人回复
-
快农历七夕节了，轻松一下，男人悄悄话，女施主请不要进来。
已经有3人回复
|
data=load('920.txt'); x=data(:,1:2); y=data(:,3); beat0=[1,1,1,1]'; [beat,r,j]=nlinfit(x,y,@myfunn,beat0) 这一段是什么意思撒 发自小木虫IOS客户端 |
2楼2017-04-09 21:10:36
3楼2017-04-12 22:16:32
4楼2017-04-15 17:14:51
 |
5楼2017-08-18 04:31:32
 |
6楼2017-08-18 19:03:06
 |
7楼2017-08-19 03:15:01
 |
8楼2017-08-19 08:55:29
|
|
9楼2017-08-20 00:02:18
 |
10楼2017-08-20 07:08:28
