伤痕累累 发表于 2011-12-30 16:38

帮我看看我的有限差分法的程序吧

我自己写得差分程序,里面有错误。高手帮我看看吧。
%%% chafenfa
function cff
clear all
clc
T=1;a=0.01;dt=0.01;m=100;h=1/m;
%Wr=100;dw=1/100;
=cfm(h,dt,a,T,m);

plot(Tt,W);

function =cfm(h,dt,a,T,m)

tol=1e-6;
N=1000;
%步长t
t=0:dt:T;%时间离散
%wr=0:dw:Wr;%角频率离散
n=length(t)
for i=2:m+2% 整体向右移2个单位,因为有-1项
    x(i)=(i-1)*h;
end
%=meshgrid(x,t);
w=zeros(n+2,m+2);
v=1.2;vf2=0.001;v12=28838;F=10;wr=2*pi;
k=1;
%while k<N %2*v/(4*h*t)=6e-5 ((v^2-1)/h^2)=4400
%边界及初始条件
for i=2:n+2
    w(2,i)=0;
    w(m+2,i)=0;
end
for j=2:m+2
    w(j,2)=a*x(j)*(1-x(j));
    w(j,n+2)=0;
end
for i=3:n+1
    for j=2:m+1
         %w(j,3)=-dt^2/2*(1.5*v12/(4*h^4)*(w(j+1,2)-2*w(j,2)+w(j-1,2))*(w(j+1,2)-w(j-1,2))^2-6e-5*(w(j+1,3)-w(j+1,1)-w(j-1,3)+w(j-1,1))+4400*(w(j+1,2)-2*w(j,2)+w(j-1,2))-F*cos(wr*dt))+w(j,2);
         w(j,1)=w(j,3);
         w(m+2,i)=w(m+1,i);         
         w(j,i)=dt^2/2*(6e-5*(w(j+1,i+1)-w(j+1,i-1)-w(j-1,i+1)+w(j-1,i-1))+4400*(w(j+1,i)-2*w(j,i)+w(j-1,i))-F*cos(wr*dt)-1.5*v12/(4*h^4)*(w(j+1,i)-2*w(j,i)+w(j-1,i))*(w(j+1,i)-w(j-1,i))^2)+1/2*(w(j,i+1)+w(j,i-1));
   
    end         
end
%k=k+1;
%end
yu=length(w);
W=w(2:m+2,2:m+2);
yu2=length(W)
Tt=t;
%if k==N+1
%      fid = fopen('FDresult.txt', 'wt');
%      fprintf(fid,'超过最大迭代次数,求解失败!');
%       fclose(fid);
%end

伤痕累累 发表于 2012-2-12 16:16

求高手指点啊

伤痕累累 发表于 2012-2-24 15:26

有限差分法的步长是应该怎么去取呢?
页: [1]
查看完整版本: 帮我看看我的有限差分法的程序吧