三次样条差值-matlab通用程序

数值分析课程的大作业，教材《数值分析》李乃成.梅立泉clearx=input('请按照格式[x1,x2,x3...]格式输入y=f(x)函数已知点的横坐标xi=');%三次样条差值函数y=input('请按照格式[y1,y2,y3...]格式输入y=f(x)函数已知点对应的纵坐标yi=');n=size(x,2);%特别注意，matlab中的矩阵编号是从1开始的，而教材上的矩阵编号是从0开始的fork=2:n%计算h(i)h(k)=x(k)-x(k-1);endfork=2:(n-1)%计算μ和λmu(k)=h(k)/(h(k)+h(k+1));lumbda(k)=1-mu(k);endfork=2:(n-1)d(k)=6*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k))/(h(k)+h(k+1));%计算diendfprintf('边界条件类型选择:\n1.已知f(a)和f(b)的二阶导数\n2.已知f(a)和f(b)的一阶导数\n3.y=f(x)是以T=b-a为周期的周期函数\n');m=input('请输入对应序号：');ifm==1M(1)=input('请输入f(a)的二阶导数值：');M(n)=input('请输入f(b)的二阶导数值：');A=zeros(n-2,n-2);fork=1:(n-3)%构造追赶法所需的A和bA(k,k)=2;A(k,k+1)=lumbda(k+1);A(k+1,k)=mu(k+2);endA(n-2,n-2)=2;b=zeros(n-2,1);fork=2:(n-3)b(k,1)=d(k+1);endb(1,1)=d(2)-mu(2)*M(1);b(n-2,1)=d(n-1)-lumbda(n-1)*M(n);N=ZhuiGanFa(A,b);%利用追赶法求解fork=2:(n-1)M(k)=N(k-1,1);endelseifm==2y0=input('请输入f(a)的一阶导数值：');yn=input('请输入f(b)的一阶导数值：');d(1)=6*((y(2)-y(1))/h(2)-y0)/h(2);d(n)=6*(yn-(y(n)-y(n-1))/h(n))/h(n);A=zeros(n,n);fork=2:(n-2)%构造追赶法所需的A和bA(k,k)=2;A(k,k+1)=lumbda(k);A(k+1,k)=mu(k+1);endA(1,1)=2;A(1,2)=1;A(2,1)=mu(2);A(n-1,n-1)=2;A(n,n)=2;A(n,n-1)=1;A(n-1,n)=lumbda(n-1);b=zeros(n,1);fork=1:nb(k,1)=d(k);endN=ZhuiGanFa(A,b);%利用追赶法求解fork=1:nM(k)=N(k,1);endelseifm==3d(n)=6*((y(2)-y(n))/h(2)-(y(n)-y(n-1))/h(n))/(h(n)+h(2));A=zeros(n-1,n-1);fork=1:(n-2)A(k,k)=2;A(k,k+1)=lumbda(k+1);A(k+1,k)=mu(k+2);endA(n-1,n-1)=2;A(1,n-1)=mu(2);A(n-1,1)=lumbda(n);b=zeros(n-1,1);fork=1:(n-1)b(k,1)=d(k+1);endN=LU_fenjieqiuxianxingfangcheng(A,b);%利用LU分解求解线性方程组fork=1:(n-1)M(k+1)=N(k,1);endM(1)=M(n);elsefprintf('您输入的序号不正确');endM;fork=1:(n-1)clearS1symsXS1=(x(k+1)-X)^3*M(k)/(6*h(k+1))+(X-x(k))^3*M(k+1)/(6*h(k+1))+...(y(k)-h(k+1)^2*M(k)/6)*(x(k+1)-X)/h(k+1)+...(y(k+1)-h(k+1)^2*M(k+1)/6)*(X-x(k))/h(k+1);fprintf('当%d=X=%d时\n',x(k),x(k+1));S=expand(S1)end----------------------------------functionx=LU_fenjieqiuxianxingfangcheng(A,b)%LU分解求解线性方程n=size(A,1);forj=1:nu(1,j)=A(1,j);endfori=2:nl(i,1)=A(i,1)/u(1,1);endfori=2:(n-1)clearSUM1SUM1=0;fork=1:(i-1)SUM1=SUM1+l(i,k)*u(k,i);endu(i,i)=A(i,i)-SUM1;forj=(i+1):nclearSUM2SUM2=0;fork=1:(i-1)SUM2=SUM2+l(i,k)*u(k,j);endu(i,j)=A(i,j)-SUM2;clearSUM3SUM3=0;fork=1:(i-1)SUM3=SUM3+l(j,k)*u(k,i);endl(j,i)=(A(j,i)-SUM3)/u(i,i);endendclearSUM4SUM4=0;fork=1:(n-1)SUM4=SUM4+l(n,k)*u(k,n);endu(n,n)=A(n,n)-SUM4;fori=1:nl(i,i)=1;endl;u;%---------------下面是利用LU分解求解线性方程组Ax=b------%y(1,1)=b(1,1);%求解Ly=bfori=2:nclearSUM5SUM5=0;forj=1:(i-1)SUM5=SUM5+l(i,j)*y(j,1);endy(i,1)=b(i,1)-SUM5;endx(n,1)=y(n,1)/u(n,n);%求解Ux=yfori=(n-1):(-1):1clearSUM5SUM5=0;forj=(i+1):nSUM5=SUM5+u(i,j)*x(j,1);endx(i,1)=(y(i,1)-SUM5)/u(i,i);end--------------------------------functionx=ZhuiGanFa(A,d)%追赶法求解线性方程组的函数n=size(A,1);u(1)=A(1,1);y(1)=d(1,1);fork=2:nclearl1l1=A(k,k-1)/u(k-1);u(k)=A(k,k)-l1*A(k-1,k);y(k)=d(k,1)-l1*y(k-1);endx(n,1)=y(n)/u(n);form=(n-1):(-1):1x(m,1)=(y(m)-A(m,m+1)*x(m+1,1))/u(m);end