下面是求解结构动力学方程的精细积分计算程序：

精细积分的优点见钟万勰教授的论著。

---

%clear;

%A=zeros(2);
%C=A;
%D=[0.5,0;0,1];
%B=[-6,2;2,-4];

%定义质量矩阵
M = [1 0 0;0 1 0;0 0 6];
%定义刚度矩阵
K = [10 1 0; 1 10 1;0 1 10];

%定义比例阻尼
alpha = 0.01;
bita = 0.02;
G = alpha*M+bita*K;

iM = inv(M);

A = -0.5*iM*G;
B = 0.25*G*iM*G-K;
C = -0.5*G*iM;
D = iM;

%
%A=zeros(3);
%C=A;
%D=[0.5 0 0;0 1 0;0 0 2];
%B=[-6 2 2;2 -4 2; 2 2 -5];

%定义初始时刻的外力向量
%f0=[0;0;0;10];
f0=[0 0 0 10 1 20]';

f1=zeros(size(f0));

%形成哈密尔顿矩阵
H=[A,D;B,C];

%计算的结束时间
tf=20; 

%下面的定义是为了验证在不同步长情况下，精细积分的精度之间的差别
step=[2,0.5,0.1]; % different step size

%这一般是精细积分计算的缺省值
N=20;

%figure;
grid;
hold on;

str=['o','x','b-'];

for jj = 1:3
   [jj tf/step(jj)]  %time,vector的维数和tf/step(jj)的值相等
   [time,vector]=jxjf1(H,f0,f1,step(jj),tf,N);
   plot(time,vector(:,1),'r-');
   plot(time,vector(:,2),'g-');
   plot(time,vector(:,3),'b-');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [time,matrix]=jxjf1(H,f0,f1,time_step,tf,N)
%
% time:计算的时间序列,一列
% matrix:计算的结果，矩阵，按照列分为两半，前一半为节点位移的响应，前一半为节点动量的响应，
%

    %定义与哈密尔顿矩阵同阶的单位矩阵
    I=eye(size(H));
   
    %需要计算增扩矩阵的逆
    iH=inv(H);
   
    %精细积分的输入为：步长、一个正整数、增扩矩阵、外力向量初值
    %PIM begin 精细积分方法开始
    dt=time_step/2^N;  %将每一个步长分为2的整幂次方份
   
    % Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12)/2;
    temp = H*dt;
    temp1 = temp*temp;
    Ta=temp+temp1*(I+temp/3+temp1/12)/2;
   
    %对每一份计算
    for iter=1:N  
        Ta=2*Ta+Ta*Ta;  %化简成Ta=（2+Ta）*Ta有些问题 huhu
    end

    %定义一个新变量
    T=I+Ta; 
   
    %vk=[0;0;0;0];  %定义一个列向量
    vk=zeros(max(size(H)),1);  %定义一个列向量 huhu
   
    number = tf/time_step;

    %number
    matrix = [];

    for iter=1:number;
       
        matrix = [matrix vk];
       
        %iter
        t(iter)=time_step*(iter-1);
        %v(iter)=vk(1);
        %vk=T*(vk+iH*(f0+iH*f1))-iH*(f0+iH*f1+f1*step(jj));
        temp = iH*(f0+iH*f1);
        vk=T*(vk+temp)-temp-iH*f1*time_step;
    end

    time = t(1:number)';
    %vector = v;

    matrix = matrix';
end

---