声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2079|回复: 2

[应用数学] 帮我看看我的有限差分法的程序吧

[复制链接]
发表于 2011-12-30 16:38 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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
回复
分享到:

使用道具 举报

 楼主| 发表于 2012-2-12 16:16 | 显示全部楼层
求高手指点啊
 楼主| 发表于 2012-2-24 15:26 | 显示全部楼层
有限差分法的步长是应该怎么去取呢?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 20:46 , Processed in 0.057662 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表