代码之家  ›  专栏  ›  技术社区  ›  user3600497

为什么我的解决方案没有腐烂

  •  0
  • user3600497  · 技术社区  · 11 年前

    我用正时-中心-空间有限差分法数值求解热方程。

    一旦空间步长足够小,我的解就会停止衰减。为什么?

    我似乎想不通,我真的很感谢你的帮助。

    这是我的代码:

    TI=0; %Intial time
    TF=1; %Final Time
    
    sigma=2;
    dx=0.01;
    dt=(dx^2)/(5*sigma); %Ensure the criteria on r is met by chooding delta t in this way
    
    r=sigma*dt/(dx^2);
    
    x=-1:dx:1;
    phi=sin(pi*x); %Initial cpndition
    
    old=phi; %To be used in the algorithm
    new=zeros(1,numel(x));
    
    timesteps=(TF-TI)/dt; %Classically called n
    timesteps=int8(timesteps); %Ensure this number is an integer so it doesnt make matlab mad
    spacesteps=numel(x);
    
    
    M=zeros(timesteps,spacesteps);
    
    M(1,:)=phi;  %Going to put all my computed time steps into a matrix.
    
    
    
    for i=2:timesteps
    
        %Now take dx space steps
    
        for j=2:spacesteps-1
    
            new(1)=0;
            new(end)=0;
            new(j)=old(j) + r*(old(j+1)-2*old(j)+old(j-1));
    
        end
    
        M(i,:)=new;
    
        old=new;
        new=zeros(1,numel(x));
    
    end
    
    
    DIM_M=size(M);
    
    [X,T]=meshgrid(linspace(-1,1,DIM_M(2)),linspace(0,TF,DIM_M(1)));
    
    figure(1)
    surf(X,T,M);
    xlabel('x')
    ylabel('t')
    title('Numerical Solution')
    shading interp
    
    AS=exp(-sigma*pi^2*T).*sin(pi*X);
    
    figure(2)
    surf(X,T,AS)
    xlabel('x')
    ylabel('t')
    title('Actual Solution')
    shading interp
    
    Error=AS-M;
    
    figure(3);
    surf(X,T,Error)
    shading interp
    
    1 回复  |  直到 11 年前
        1
  •  1
  •   Ben Voigt    11 年前

    选择导致错误行为的几行。

    TI=0;
    TF=1;
    
    sigma=2;
    dx=0.01;
    dt=(dx^2)/(5*sigma); // .01^2 == .0001, 5*sigma == 10, dt := .00001
    
    timesteps=(TF-TI)/dt; // TF-TI == 1, 1/.00001 == 10000, timesteps := 10000
    timesteps=int8(timesteps);
    

    最后一行出现超出范围的情况。MATLAB剪辑在类型转换中溢出到最接近的可表示值,因此实际上 timesteps = 127 和你不变的 TF == TI + dt * timesteps; 被违反。模拟将远远低于 TF .