数值分析课程的大作业，教材《数值分析》李乃成.梅立泉

clear

x=input(‘请按照格式[x1,x2,x3…]格式输入y=f(x)函数已知点的横座标xi=’);   %三次样条差值函数

y=input(‘请按照格式[y1,y2,y3…]格式输入y=f(x)函数已知点对应的纵座标yi=’);

n=size(x,2);      %特别注意，matlab中的矩阵编号是从1开始的，而教材上的矩阵编号是从0开始的

for k=2:n           %计算h(i)

h(k)=x(k)-x(k-1);

end

for k=2:(n-1)                    %计算μ和λ

mu(k)=h(k)/(h(k)+h(k+1));

lumbda(k)=1-mu(k);

end

for k=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));            %计算di

end

fprintf(‘边界条件类型选择:\n1.已知f(a)和f(b)的二阶导数\n2.已知f(a)和f(b)的一阶导数\n3.y=f(x)是以T=b-a为周期的周期函数\n’);

m=input(‘请输入对应序号：’);

if m==1

M(1)=input(‘请输入f(a)的二阶导数值：’);

M(n)=input(‘请输入f(b)的二阶导数值：’);

A=zeros(n-2,n-2);

for k=1:(n-3)             %构造追赶法所需的A和b

A(k,k)=2;

A(k,k+1)=lumbda(k+1);

A(k+1,k)=mu(k+2);

end

A(n-2,n-2)=2;

b=zeros(n-2,1);

for k=2:(n-3)

b(k,1)=d(k+1);

end

b(1,1)=d(2)-mu(2)*M(1);

b(n-2,1)=d(n-1)-lumbda(n-1)*M(n);

N=ZhuiGanFa(A,b);            %利用追赶法求解

for k=2:(n-1)

M(k)=N(k-1,1);

end

elseif m==2

y0=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);

for k=2:(n-2)                      %构造追赶法所需的A和b

A(k,k)=2;

A(k,k+1)=lumbda(k);

A(k+1,k)=mu(k+1);

end

A(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);

for k=1:n

b(k,1)=d(k);

end

N=ZhuiGanFa(A,b);            %利用追赶法求解

for k=1:n

M(k)=N(k,1);

end

elseif m==3

d(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);

for k=1:(n-2)

A(k,k)=2;

A(k,k+1)=lumbda(k+1);

A(k+1,k)=mu(k+2);

end

A(n-1,n-1)=2;

A(1,n-1)=mu(2);

A(n-1,1)=lumbda(n);

b=zeros(n-1,1);

for k=1:(n-1)

b(k,1)=d(k+1);

end

N=LU_fenjieqiuxianxingfangcheng(A,b);       %利用LU分解求解线性方程组

for k=1:(n-1)

M(k+1)=N(k,1);

end

M(1)=M(n);

else

fprintf(‘您输入的序号不正确’);

end

M;

for k=1:(n-1)

clear S1

syms X

S1=(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=

S=expand(S1)

end

———————————-

function x=LU_fenjieqiuxianxingfangcheng(A,b)      %LU分解求解线性方程

n=size(A,1);

for j=1:n

u(1,j)=A(1,j);

end

for i=2:n

l(i,1)=A(i,1)/u(1,1);

end

for i=2:(n-1)

clear SUM1

SUM1=0;

for k=1:(i-1)

SUM1=SUM1+l(i,k)*u(k,i);

end

u(i,i)=A(i,i)-SUM1;

for j=(i+1):n

clear SUM2

SUM2=0;

for k=1:(i-1)

SUM2=SUM2+l(i,k)*u(k,j);

end

u(i,j)=A(i,j)-SUM2;

clear SUM3

SUM3=0;

for k=1:(i-1)

SUM3=SUM3+l(j,k)*u(k,i);

end

l(j,i)=(A(j,i)-SUM3)/u(i,i);

end

end

clear SUM4

SUM4=0;

for k=1:(n-1)

SUM4=SUM4+l(n,k)*u(k,n);

end

u(n,n)=A(n,n)-SUM4;

for i=1:n

l(i,i)=1;

end

l;

u;

%—————下面是利用LU分解求解线性方程组Ax=b——%

y(1,1)=b(1,1);%求解Ly=b

for i=2:n

clear SUM5

SUM5=0;

for j=1:(i-1)

SUM5=SUM5+l(i,j)*y(j,1);

end

y(i,1)=b(i,1)-SUM5;

end

x(n,1)=y(n,1)/u(n,n);%求解Ux=y

for i=(n-1):(-1):1

clear SUM5

SUM5=0;

for j=(i+1):n

SUM5=SUM5+u(i,j)*x(j,1);

end

x(i,1)=(y(i,1)-SUM5)/u(i,i);

end

——————————–

function x=ZhuiGanFa(A,d)   %追赶法求解线性方程组的函数

n=size(A,1);

u(1)=A(1,1);

y(1)=d(1,1);

for k=2:n

clear l1

l1=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);

end

x(n,1)=y(n)/u(n);

for m=(n-1):(-1):1

x(m,1)=(y(m)-A(m,m+1)*x(m+1,1))/u(m);

end