求助姜启源<数学建模案例选集>中的一个程序

大白菜

程序原文:
global wij;
fid=fopen('data.txt','w');
T=0.1;t0=0;tn=0;
dx=0.01;dy=0.01;dt=0.001;c=2000;
nx=60;ny=60;
u=zeros(nx,ny);un=zeros(nx,ny);
vn=zeros(nx,ny);v=zeros(nx,ny);
u(nx/2,ny/2)=1e-4;
u(nx/2-1,ny/2)=1e-4/2;u(nx/2+1,ny/2)=1e-4/2;
u(nx/2,ny/2-1)=1e-4/2;u(nx/2,ny/2+1)=1e-4/2;
u(nx/2-1,ny/2-1)=1e-4/2;u(nx/2-1,ny/2+1)=1e-4/2;
u(nx/2+1,ny/2-1)=1e-4/2;u(nx/2+1,ny/2+1)=1e-4/2;
%
while(tn<=T)
    t0=tn;tn=tn+dt;ts=[t0,tn];
      for i=2:(nx-1)
          for j=2:(ny-1)
              wij=c*c*(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx*dx)+c*c*(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy*dy);
              w0=[u(i,j),v(i,j)];
              [t,w]=ode23('wave',ts,w0);n1=length(t);
              un(i,j)=w(n1,1);vn(i,j)=w(n1,2);
          end
      end
      u=un;v=vn;
%
       for i=1:nx
           for j=1:ny
               fprintf(fid,'%.3e\n',u(i,j));
           end
       end
%
end
surfl(un)
shading interp
colormap(gray);
fclose(fid);
=====================================================================================
我的结果是 jieguo.fig (68.53 KB, 下载次数: 10)

2008-8-14 09:32 上传

点击文件名下载附件

显然不对,请问是哪里出问题了?

大白菜

谢谢messenger 对上个贴子的协助

messenger

贴个图片吧，fig图打开什么也看不到呀

大白菜

贴图我的不会啊

sigma665

把图片保存为jpg格式

大白菜

贴图如下,请看看

[ 本帖最后由 大白菜 于 2008-8-14 15:20 编辑 ]

大白菜

好的,谢谢斑竹:

大白菜

这个结果差太远了,我怀疑是un的值不对，但是无法验证

messenger

你这个图片上什么也没有，和fig上看到的情况一样，但点fig上好象能点到东西，就好象图像被画成了白色一样，但是我改图像颜色也没显示出东西来，是不是你画图有问题呀

大白菜

原帖由 messenger 于 2008-8-14 15:22 发表

                               
登录/注册后可看大图

你这个图片上什么也没有，和fig上看到的情况一样，但点fig上好象能点到东西，就好象图像被画成了白色一样，但是我改图像颜色也没显示出东西来，是不是你画图有问题呀

==================================================================
我没有觉得啊?你能提示一下么?

sigma665

图上除了坐标，没有实质内容

messenger

就是点击一下图的中央，可以看到有几个轮廓点，就象点到其它的三维图一样

也可能像版主说的，什么也没有，你还是自己检查一下吧。

大白菜

谢谢,问题是我那是二维图像,怎么会画出三维图啊?

messenger

本来就是三维图像，有三个坐标轴呀

不过，看起来像二维的图像，原书上的图像也是这样，看起来像二维的，反正如果你严格按原书来做肯定没错。

大白菜

谢谢你,我再查查抄写的时候有没有错,我发现运行后un的值大多为NaN,这又是怎么回事?怎么都是不定值???
====================================================================
原文和我的程序有两个差别,一个就是c我赋值了,第二就是:
wij=c*c*(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx*dx)+c*c*(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy*dy);
原文为:
wij=c*c*(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx*dx)+...c*c*(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy*dy);

=======================================================================这是怎么回事？难道是我设置的数据类型与运算结果有冲突？
请各位高手为我指点迷津 ，谢谢