浙江大学计算方法大作业第九题

结果分析：
结果较为精确，但是由于追赶法解线性方程组的过程中，由于追赶法本身的特点，会是 M0
的值产生误差，从而影响结果，当然这个误差在可接受的范围内。
源程序：
#include<stdio.h> int main(void) { FILE *fp; int i,j; double x[19]={ 0.52, 3.1, 8.0, 17.95, 28.65, 39.62, 50.65, 78, 104.6, 156.6, 208.6, 260.7, 312.5, 364.4, 416.3, 468, 494, 507, 520}; double y[19]={5.28794, 9.4, 13.84, 20.2, 24.9, 28.44, 31.1, 35, 36.5, 36.6, 34.6, 31.0, 26.34, 20.9, 14.8, 7.8, 3.7, 1.5, 0.2}; double a[12]={2, 4, 6, 12, 16, 30, 60, 110, 180, 280, 400, 515}; double b[12]; double ya=1.86548, yb= -0.046115; double h[19],u[19],l[19],g[19]; double bt[19],af[19],f[19],m[19]; for(i=1;i<=18;i++) h[i]=x[i]-x[i-1]; for(i=1;i<=17;i++){ u[i]=h[i]/(h[i]+h[i+1]); l[i]=1-u[i]; g[i]=6/(h[i]+h[i+1])*((y[i+1]-y[i])/h[i+1]-(y[i]-y[i-1])/h[i]); } g[0]=6/h[1]*((y[1]-y[0])/h[1]-ya); g[18]=6/h[18]*(yb-(y[18]-y[17])/h[18]); u[18]=1.0; l[0]=1.0; af[0]=2.0; for(i=1;i<=18;i++){ bt[i-1]=l[i-1]/af[i-1]; af[i]=2.0-u[i]*bt[i-1]; } f[0]=g[0]/2.0; for(i=1;i<=18;i++) f[i]=(g[i]-u[i]*f[i-1])/(2.0-u[i]*bt[i-1]);
( x [ xi 1 , xi ], i 1, 2, ……，n) 对于求M i，则需要用到追赶法，详见课本P44-P45
计算结果：
s(2)=7.825123 s(4)=10.481311 s(6)=12.363477 s(12)=16.575574 s(16)=19.091594 s(30)=25.386402 s(60)=32.804283 sБайду номын сангаас110)=36.647878 s(180)=35.917148 s(280)=29.368462 s(400)=16.799784 s(515)=0.542713
题目：9：已知直升飞机旋转及以外形曲线轮廓线上的某些型值点（见表）及端点处的一
阶导数值试计算该曲线上横坐标为 2,4,6,12,16,30,60,110,180,280,400,515 处点的纵坐标（要求该曲线具有二阶光滑度） 。 k xk yk k xk yk k xk yk 解： 0 0.52 5.28794 7 78 35 14 416.3 14.8 1 3.1 9.4 8 104.6 36.5 15 468 7.8 2 8.0 13.84 9 156.6 36.6 16 494 3.7 3 17.95 20.2 10 208.6 34.6 17 507 1.5 4 28.65 24.9 11 260.7 31.0 18 520 0.2 5 39.62 28.44 12 312.5 26.34 6 50.65 31.1 13 364.4 20.9
m[18]=f[18]; for(i=17;i>=0;i--) m[i]=f[i]-bt[i]*m[i+1]; for(j=0;j<=11;j++){ for(i=1;i<=18;i++){ if(a[j]>=x[i-1]&&a[j]<=x[i]){ b[j]=m[i-1]*(x[i]-a[j])*(x[i]-a[j])*(x[i]-a[j])/(6.0*h[i]) +m[i]*(a[j]-x[i-1])*(a[j]-x[i-1])*(a[j]-x[i-1])/(6.0*h[i]) +(y[i-1]-m[i-1]*h[i]*h[i]/6.0)*(x[i]-a[j])/h[i] +(y[i]-m[i]*h[i]*h[i]/6.0)*(a[j]-x[i-1])/h[i]; break; } else continue; } } fp=fopen("9.txt","w"); for(j=0;j<=11;j++) fprintf(fp,"s(%.0f)=%f\n",a[j],b[j]); fclose(fp); }
思路步骤：
利用三次样条插值法：
按照条件整理后用M i 表示si ( x)的公式 有si ( x) M i 1 ( xi x)3 ( x xi 1 )3 M x x M x xi 1 Mi ( yi 1 i 1 hi2 ) i ( yi i hi2 ) 6hi 6hi 6 hi 6 hi