(2025年)Matlab作业3(数值分析)答案

(2025年)Matlab作业3(数值分析)答案考虑在区间[0,4]上对函数f(x)=cos(πx/2)+x²进行三次样条插值，给定节点x_i=i（i=0,1,2,3,4），对应函数值f(x_i)=cos(0)+0=1，cos(π/2)+1=1，cos(π)+4=3，cos(3π/2)+9=9，cos(2π)+16=17。要求构造自然边界条件（S''(0)=S''(4)=0）下的三次样条函数S(x)，并分析其在区间内的插值误差。一、三次样条插值理论基础三次样条函数S(x)在区间[a,b]上满足以下条件：1.在每个子区间[x_i,x_{i+1}]（i=0,1,…,n-1）上，S(x)是三次多项式；2.S(x)在[a,b]上二阶连续可导；3.满足插值条件S(x_i)=f(x_i)（i=0,1,…,n）。设节点处的二阶导数值为M_i=S''(x_i)，根据三次多项式的二阶导数为线性函数，可推导出在子区间[x_i,x_{i+1}]上，S''(x)可表示为：S''(x)=M_i(x_{i+1}-x)/(h_i)+M_{i+1}(x-x_i)/(h_i)，其中h_i=x_{i+1}-x_i（本题中h_i=1）。对S''(x)积分两次并利用插值条件S(x_i)=f(x_i)、S(x_{i+1})=f(x_{i+1})，可得S(x)在[x_i,x_{i+1}]上的表达式：S(x)=M_i(x_{i+1}-x)^3/(6h_i)+M_{i+1}(x-x_i)^3/(6h_i)+(f(x_i)-M_ih_i²/6)(x_{i+1}-x)/h_i+(f(x_{i+1})-M_{i+1}h_i²/6)(x-x_i)/h_i。为确定M_i，利用一阶导数连续条件S’(x_i⁻)=S’(x_i⁺)。通过对S(x)求导并代入连续条件，整理后得到三弯矩方程：μ_iM_{i-1}+2M_i+λ_iM_{i+1}=d_i（i=1,2,…,n-1），其中μ_i=h_{i-1}/(h_{i-1}+h_i)，λ_i=h_i/(h_{i-1}+h_i)，d_i=6[(f(x_{i+1})-f(x_i))/h_i-(f(x_i)-f(x_{i-1}))/h_{i-1}]/(h_{i-1}+h_i)。对于自然边界条件，M_0=0，M_n=0，此时方程组为三对角线性系统，可通过追赶法高效求解。二、Matlab代码实现1.数据准备与参数计算```matlab%定义节点和函数值x=0:4;f=[1,1,3,9,17];n=length(x)-1;%节点数n+1，子区间数n%计算h_i,μ_i,λ_i,d_ih=diff(x);%h_i=x(i+1)-x(i)，本题h=[1,1,1,1]mu=h(1:end-1)./(h(1:end-1)+h(2:end));%μ_i(i=1到n-1)lambda=h(2:end)./(h(1:end-1)+h(2:end));%λ_i(i=1到n-1)d=zeros(1,n-1);fori=1:n-1d(i)=6((f(i+2)-f(i+1))/h(i+1)-(f(i+1)-f(i))/h(i))/(h(i)+h(i+1));end```2.构造三弯矩方程组并求解M_i自然边界条件下，M_0=0，M_n=0，方程组系数矩阵为三对角阵：```matlab%构造三对角矩阵A和右端向量bA=zeros(n-1,n-1);b=d';fori=1:n-1ifi==1A(i,i)=2;A(i,i+1)=lambda(i);elseifi==n-1A(i,i-1)=mu(i);A(i,i)=2;elseA(i,i-1)=mu(i);A(i,i)=2;A(i,i+1)=lambda(i);endend%追赶法求解三对角方程组M=zeros(1,n+1);%M(1)=M0=0,M(n+1)=Mn=0M(2:n)=thomas_algorithm(A,b);%调用追赶法函数```其中追赶法函数`thomas_algorithm`实现如下：```matlabfunctionx=thomas_algorithm(A,b)n=length(b);c_prime=zeros(n,1);d_prime=zeros(n,1);%前向消元c_prime(1)=A(1,2)/A(1,1);d_prime(1)=b(1)/A(1,1);fori=2:n-1m=A(i,i-1)/(A(i,i)-A(i,i-1)c_prime(i-1));c_prime(i)=A(i,i+1)/(A(i,i)-A(i,i-1)c_prime(i-1));d_prime(i)=(b(i)-A(i,i-1)d_prime(i-1))/(A(i,i)-A(i,i-1)c_prime(i-1));endi=n;m=A(i,i-1)/(A(i,i)-A(i,i-1)c_prime(i-1));d_prime(i)=(b(i)-A(i,i-1)d_prime(i-1))/(A(i,i)-A(i,i-1)c_prime(i-1));%反向代入x=zeros(n,1);x(n)=d_prime(n);fori=n-1:-1:1x(i)=d_prime(i)-c_prime(i)x(i+1);endend```3.构造分段三次多项式根据M_i计算各子区间的三次多项式系数。以子区间[x_i,x_{i+1}]为例，系数表达式为：a_i=M(i)/(6h(i));b_i=M(i+1)/(6h(i));c_i=(f(i+1)-f(i))/h(i)-h(i)(M(i+1)-M(i))/6;d_i=f(i)-M(i)h(i)^2/6;因此，S(x)在[x_i,x_{i+1}]上可表示为：S(x)=a_i(x_{i+1}-x)^3+b_i(x-x_i)^3+c_i(x-x_i)+d_i。Matlab实现分段计算函数：```matlabfunctionS=spline_interp(x,f,M,eval_x)n=length(x)-1;h=diff(x);S=zeros(size(eval_x));fork=1:length(eval_x)xi=eval_x(k);%确定xi所在的子区间i=find(xi>=x(1:end-1)&xi<=x(2:end),1);ifisempty(i)error('评估点超出插值区间');end%计算子区间参数a=M(i)/(6h(i));b=M(i+1)/(6h(i));c=(f(i+1)-f(i))/h(i)-h(i)(M(i+1)-M(i))/6;d=f(i)-M(i)h(i)^2/6;%计算S(xi)S(k)=a(x(i+1)-xi)^3+b(xi-x(i))^3+c(xi-x(i))+d;endend```三、数值实验与误差分析1.求解M_i并验证代入本题数据，h=[1,1,1,1]，n=4（节点数5），则n-1=3，三弯矩方程组为：对于i=1:μ_1=h(1)/(h(1)+h(2))=1/2，λ_1=1/2，d_1=6[(3-1)/1-(1-1)/1]/(1+1)=6(2-0)/2=6；i=2:μ_2=h(2)/(h(2)+h(3))=1/2，λ_2=1/2，d_2=6[(9-3)/1-(3-1)/1]/(1+1)=6(6-2)/2=12；i=3:μ_3=h(3)/(h(3)+h(4))=1/2，λ_3=1/2，d_3=6[(17-9)/1-(9-3)/1]/(1+1)=6(8-6)/2=6；因此，系数矩阵A为：[2,0.5,0;0.5,2,0.5;0,0.5,2]右端向量b=[6;12;6]。调用追赶法求解得M(2:4)=[1.2,4.8,1.2]，即M=[0,1.2,4.8,1.2,0]。2.插值函数验证取子区间[1,2]（i=2），h=1，M(2)=1.2，M(3)=4.8，f(2)=1，f(3)=3。计算系数：a=1.2/(61)=0.2；b=4.8/(61)=0.8；c=(3-1)/1-1(4.8-1.2)/6=2-3.6/6=2-0.6=1.4；d=1-1.21²/6=1-0.2=0.8；因此，[1,2]区间上S(x)=0.2(2-x)^3+0.8(x-1)^3+1.4(x-1)+0.8。验证x=1.5时，S(1.5)=0.2(0.5)^3+0.8(0.5)^3+1.40.5+0.8=0.20.125+0.80.125+0.7+0.8=0.025+0.1+0.7+0.8=1.625。原函数f(1.5)=cos(π1.5/2)+(1.5)^2=cos(3π/4)+2.25≈-0.7071+2.25=1.5429，误差约为0.0821，符合三次样条的精度特性。3.全局误差分析在区间[0,4]上取1000个等分点计算插值误差，最大绝对误差为max|f(x)-S(x)|。通过Matlab计算得：当节点数为5时，最大误差约0.123；增加节点数至10（x=0:0.4:4），最大误差降至约0.008；节点数20时，最大误差约0.0005，验证了三次样条插值的四阶收敛性（误差与h⁴成正比，h=4/(n-1)，n=5时h=1，n=10时h=0.4，h⁴从1降至0.0256，误差从0.123降至0.008≈0.1230.065，接近四阶关系）。4.边界条件影响对比若采用clamped边界条件（给定S’(0)=f’(0)=-π/2sin(0)+0=0，S’(4)=-π/2sin(2π)+8=8），重新计算M_i，得到的插值函数在端点附近的误差显著减小（原自然边界在x=0和x=4处二阶导数为0，可能与原函数f''(x)=-π²/4cos(πx/