关于中心差分法请教各位大侠指教
本帖最后由 xzning123 于 2012-3-31 09:46 编辑这个是我写的程序，望大家指教
以下为程序代码，%%%%central difference methrod
clc
clear all
close all
%%%%intial values of the system
m=0.2533;
c=0.1592;
k=10;
%%%%the expression of excitation
syms t
p=10*sin(pi*t/0.6);
%%%%intial caculation
t0=0.6; %the end of the excitation
total=1.0; %totla caculation time
dt=0.1; %time step
n=total/dt;
n1=t0/dt;
u1=0;
v1=0;
a1=(subs(p,0)-c*v1-k*u1)/m;
uu=u1-dt*v1+0.5*dt^2*a1;
keq=m/dt^2+c/2/dt;
x=m/dt^2-c/2/dt;
y=k-2*m/dt^2;
%%%%form the vector of the%%%%
u(1)=uu;
u(2)=u1;
v(1)=v1;
a(1)=a1;
%%%%for every time step%%%%
for i=2:1:n+1
if i<n1+1
pt=p;
else
pt=0;
end
peq(i)=subs(pt,(i-1)*dt)-x*u(i-1)-y*u(i);
u(i+1)=peq(i)/keq;
end
%%%%to figure out the caculation result
t1=0:dt:total;
plot(t1',u(:,2:n+2)) %the value corresponding to u(1) is not real
此程序参考的是《结构动力学：理论及其在地震工程中的应用》一书的中心差分法写的，其间遇到的问题是：
1、力的表达式是用符号函数写的，这样运行比较慢，是否有其他写法；
2、由于力是分段的，这样调用力的时候就用到了if语句，降低了运算效率，有没有其他写法
本人是初学者，诚请各位大侠帮助，小弟不胜感激
页:
[1]