非线性振动 matlab,非线性振动肿么编写matlab程序

function [y0,C,x]=Spline_interp_1(X,Y,s0,sN,x0)

%X,Y是已知插值点坐标

%s0,sN是两端点的一次导数值

%x0是插值点

%y0是三次样条函数在x0处的值

%C是分段三次样条函数的系数,按升幂排列,区间[Xj,Xj+1]上的表达式为C(j,1)+C(j,2)*(x-Xj)++C(j,3)*(x-Xj)^2+C(j,4)*(x-Xj)^3.

%x是用追赶法依下标顺利求得的插值节点X,Y处的二阶导数值

N=length(X);

C=zeros(4,N-1); h=zeros(1,N-1);

mu=zeros(1,N-1); lmt=zeros(1,N-1);

d=zeros(1,N); %d表示右端函数值

h=X(1,2:N)-X(1,1:N-1);

mu(1,N-1)=1; lmt(1,1)=1;

mu(1,1:N-2)=h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1));

lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1:N-2)+h(1,2:N-1));

d(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);

d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);

d(1,2:N-1)=6*((Y(1,3:N)-Y(1,2:N-1))./h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2))./h(1,1:N-2))./(h(1,1:N-2)+h(1,2:N-1));

%追赶法解三对角方程组

bit