马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我自己写得差分程序,里面有错误。高手帮我看看吧。
%%% chafenfa
function cff
clear all
clc
T=1;a=0.01;dt=0.01;m=100;h=1/m;
%Wr=100;dw=1/100;
[W,Tt]=cfm(h,dt,a,T,m);
plot(Tt,W);
function [W,Tt]=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
%[X,Y]=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 |