微分求积法(Differential Quadrature Method)
微分求积法,简称 DQ,是一种求解常微分和偏微分方程(组)的有效方法。将常微分和偏微分方程(组)经过DQ离散成线性或非线性方程组后可通过Newmark 和 Newton-Rapuson 求解,也可以直接使用Runge-Kutta法。但是本人在使用这些方法是遇到了困难,主要的问题是边界条件的处理,之后编程和算法上,因为它的离散过程是简单的!很希望和正在或曾经使用此方法的朋友交流,并虚心请教! 参考下面的程序吧% Determination of the coefficients in DQM
% 19/03/2006 By Dr. Yang J., Modified by Mr. Xiang H.J.
% Department of BC, City University of Hong Kong, Hong Kong
% =DQC(N,M,KOP)
% KOP, Keyoption. 1--uniform space pattern; 2,3--cosine space pattern
% N--Total grid point; M--Order
function =DQC(N,M,KOP)
%1. Generation of grid point coordinates in X-axis
for i=1:N
X1(i)=(i-1)/(N-1);
X2(i)=0.5*(1-cos(pi*(i-1)/(N-1)));
end
X3(1)=0;X3(2)=0.00000001;%0.00001;
for i=3:N
X3(i)=0.5*(1-cos(pi*(i-1)/(N-1)));
end
X4(1)=0;X4(2)=0.00001;X4(N-1)=0.99999;X4(N)=1.0;
for i=3:N-2
X4(i)=0.5*(1-cos(pi*(i-2)/(N-3)));
end
X5(1)=0;X5(2)=0.00001;X5(3)=0.001;
X5(N-2)=0.999;X5(N-1)=0.99999;X5(N)=1.0;
for i=4:N-3
X5(i)=0.5*(1-cos(pi*(i-3)/(N-5)));
end
if KOP==1
X=X1;
elseif KOP==2
X=X2;
elseif KOP==3
X=X3;
elseif KOP==4
X=X4;
else
X=X5;
end
%2. Producing the weighting coefficients
%2.1 Producing the first-order weighting coefficients
C=zeros(N,N,M);
for i=1:N
Q(i)=1.0;
for j=1:N
if i~=j
Q(i)=Q(i)*(X(i)-X(j));
end
end
end
for i=1:N
for j=1:N
if i~=j
C(i,j,1)=Q(i)/(X(i)-X(j))/Q(j);
C(i,i,1)=C(i,i,1)-C(i,j,1);
end
end
end
%2.2 Producing the High-order weighting coefficients
if M>1
for k=2:M
for i=1:N
for j=1:N
if i~=j
C(i,j,k)=k*(C(i,i,k-1)*C(i,j,1)-C(i,j,k-1)/(X(i)-X(j)));
C(i,i,k)=C(i,i,k)-C(i,j,k);
end
end
end
end
end
来自:http://hi.baidu.com/hjxiang/blog/item/91e66be79543902fb93820fb.html
回复 #2 风花雪月 的帖子
谢谢!这个我已经做过了,而且已经用C重新写了一遍。现在的问题就是求解二阶微分方程了,我要用4阶龙格库塔法!正在思考。其实我倒自己写了一个,但是不是很好用。算起来很慢,不知道为什么!楼主,能不能帮着找一找更好的程序!回复 #3 su200703 的帖子
我想问一下,为什么不直接用ode45这些,命令直接求解,这个方法精度更高吗? 原帖由 su200703 于 2007-10-26 16:55 发表 http://www.chinavib.com/forum/images/common/back.gif谢谢!这个我已经做过了,而且已经用C重新写了一遍。现在的问题就是求解二阶微分方程了,我要用4阶龙格库塔法!正在思考。其实我倒自己写了一个,但是不是很好用。算起来很慢,不知道为什么!楼主,能不能帮着找 ...
rk法本来的效率就不高,龙格库塔法的程序网络上有很多,你可以借鉴,本站搜索一下就能找到 原帖由 无水1324 于 2007-10-26 23:33 发表 http://www.chinavib.com/forum/images/common/back.gif
我想问一下,为什么不直接用ode45这些,命令直接求解,这个方法精度更高吗?
ode45是变步长的,说不定楼主想用定步长或者用自己学习编程也未可知 还与人看这个帖子吗?程序看不懂,谁给发个结合实例的?
页:
[1]