解微分方程求求求助!!!!!

大白菜

function wdot=wave(t,w)
global wij;
wdot=[0.1,wij];

global wij;
fid=fopen('data.txt','w');
T=0.15;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);
============================================================================
报错为
??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.
Error in ==> D:\MATLAB6p5\toolbox\matlab\funfun\ode23.m
On line 154  ==> [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, ...

我只知道是迭代步数超过限制,那么下面那个错是什么???谢谢各位老手帮我看看,这是教材上的程序原文,该不得错哒.

messenger

wave函数你没有给全吧，我这的错误和你的错误不一样

function wdot=wave(t,w)
global wij;
wdot=[0.1,wij];

大白菜

首先谢谢,但是我还是没有搞懂.请问你的错误是什么??
再次感激

sigma665

你的微分方程是什么

大白菜

===================================================
function wdot=wave(t,w)
global wij;
wdot=[0.1,wij];

global wij;
=======================================================
我这个是从姜启源的<数学建模选集>里看到这个程序,模拟二维地震波场的,但是程序却无法运行.希望大家看看是怎么回事.

messenger

wdot=[0.1,wij]; 应该是 wdot=[0.1; wij]; 吧

大白菜

不对,主要是
===================================================
function wdot=wave(t,w)
global wij;
wdot=[0.1,wij];

global wij;
=======================================================
这砣全局变量我不知道是什么意思

messenger

查了一下原书，是wdot=[w(2), wij]'; 用wdot=[w(2), wij]'; 运行了一下没有错误。

不知道你为什么改成了wdot=[0.1,wij]; 不过就算你用wdot=[0.1; wij];也应该能算出来。

另外，global wij;是全局变量的意思。

从你的回贴来看，你把第二个global wij;也放到 wave函数里了。

第二个global wij;应该和主程序放在一起的，你可能是这里错了。

原帖由 大白菜 于 2008-8-13 15:34 发表

                               
登录/注册后可看大图

不对,主要是
===================================================
function wdot=wave(t,w)
global wij;
wdot=[0.1,wij];

global wij;
=======================================================
这砣全局 ...

[ 本帖最后由 messenger 于 2008-8-13 16:02 编辑 ]

大白菜

谢谢老哥,我改了一下,原文c没有赋值,我加上了.但是结果不对呢

messenger

这个注意到了，不过c应该只影响衰减速率，不影响图像的形状

大白菜

但是结果图像不对啊,

大白菜

你的结果是对的么?