东南大学数值分析上机实验题(下)

东南大学数值分析上机实验题(下)

数值分析上机报告

姓名： 学号：

2013年12月22日

东南大学数值分析上机实验题(下)

第四章

38.（上机题）3次样条插值函数

（1）编制求第一型3次样条插值函数的通用程序;

''

端点条件为y0=0.8，y10=0.2。用所编制程序求车门的3次样条插值函数S(x),并打印出

S(i+0.5)(i=0,1,…9)。

解：

（1）

#include <iostream.h> #include <math.h>

float x1[100],f1[100],f2[99],f3[98],m[100],a[100][101],x,d[100]; float c[99],e[99],h[99],u[99],w[99],y_0,y_n,arr,s; int i,j,k,n,q;

void selectprint(float y) {

if ((y>0)&&(y!=1)) cout<<"+"<<y; else if (y==1) cout<<"+"; else if (y<0) cout<<y; }

void printY(float y){ if (y!=0) cout<<y; }

float calculation(float x){ for (j=1;j<=n;j++) if (x<=x1[j]) {

s=(float)(f1[j-1]+c[j-1]*(x-x1[j-1])+m[j-1]/2.0*(x-x1[j-1])*(x-x1[j-1])+e[j-1]*(x-x1[j-1])*(x-x1[j-1])*(x-x1[j-1])); break; }

return s; }

void main() {

东南大学数值分析上机实验题(下)

do{

cout<<"请输入n值:"; cin>>n;

if ((n>100)||(n<1)) cout<<"请重新输入整数(1..100):"<<endl; }

while ((n>100)||(n<1));

cout<<"请输入Xi(i=0,1,...,"<<n<<"):"; for (i=0;i<=n;i++) cin>>x1[i]; cout<<endl;

cout<<"请输入Yi(i=0,1,...,"<<n<<"n):"; for (i=0;i<=n;i++) cin>>f1[i]; cout<<endl;

cout<<"请输入第一型边界条件Y'0，Y'n:"; cin>>y_0>>y_n; cout<<endl;

for (i=0;i<n;i++) h[i]=x1[i+1]-x1[i];

for (i=1;i<n;i++) u[i]=(float) (h[i-1]/(h[i-1]+h[i])); for (i=1;i<n;i++) w[i]=(float) (1.0-u[i]);

for (i=0;i<n;i++) f2[i]=(f1[i+1]-f1[i])/h[i]; //一阶差商 for (i=0;i<n-1;i++) f3[i]=(f2[i+1]-f2[i])/(x1[i+2]-x1[i]); //二阶差商 for (i=1;i<n;i++) d[i]=6*f3[i-1]; //求出所有的d[i] d[0]=6*(f2[0]-y_0)/h[0]; d[n]=6*(y_n-f2[n-1])/h[n-1]; for (i=0;i<=n;i++) for (j=0;j<=n;j++)

if (i==j) a[i][j]=2; else a[i][j]=0; a[0][1]=1; a[n][n-1]=1;

for (i=1;i<n;i++) {

a[i][i-1]=u[i]; a[i][i+1]=w[i]; }

for (i=0;i<=n;i++) a[i][n+1]=d[i];

for (i=1;i<=n;i++) //用追赶法解方程,得M[i] {

arr=a[i][i-1];

for (j=0;j<=n+1;j++)

a[i][j]=a[i][j]-a[i-1][j]*arr/a[i-1][i-1]; }

m[n]=a[n][n+1]/a[n][n];

for (i=n-1;i>=0;i--) m[i]=(a[i][n+1]-a[i][i+1]*m[i+1])/a[i][i];

东南大学数值分析上机实验题(下)

for (i=0;i<n;i++) //计算S（x）的表达式

c[i]=(float) (f2[i]-(1.0/3.0*m[i]+1.0/6.0*m[i+1])*h[i]); for (i=0;i<n;i++)

e[i]=(m[i+1]-m[i])/(6*h[i]); for (i=0;i<n;i++) {

cout<<"X属于区间["<<x1[i]<<","<<x1[i+1]<<"]时"<<endl<<endl; cout<<"S(x)="; printY(f1[i]); if (c[i]!=0) {

selectprint(c[i]); cout<<"(x"; printY(-x1[i]); cout<<")"; }

if (m[i]!=0) {

selectprint((float)(m[i]/2.0)); for (q=0;q<2;q++) {

cout<<"(x"; printY(-x1[i]); cout<<")"; } }

if (e[i]!=0) {

selectprint(e[i]); for (q=0;q<3;q++) {

cout<<"(x"; printY(-x1[i]); cout<<")"; } }

cout<<endl<<endl; }

cout<<"针对本题，计算S(i+0.5)(i=0,1,...,9):"<<endl; for (i=0;i<10;i++) {

if ((i+0.5<=x1[n])&&(i+0.5>=x1[0])) {

东南大学数值分析上机实验题(下)

calculation((float)(i+0.5));

cout<<"S("<<i+0.5<<")="<<s<<endl; }

else cout<<i+0.5<<"超出定义域"<<endl; }; cout<<endl; } （2）

编制的程序求车门的3次样条插值函数S(x)：

x属于区间[0,1] 时；S(x)=2.51+0.8(x)-0.0014861(x)(x)-0.00851395(x)(x)(x)

x属于区间[1,2] 时；S(x)=3.3+0.771486(x-1)-0.027028(x-1)(x-1)-0.00445799(x-1)(x-1)(x-1) x属于区间[2,3] 时；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)(x-2)-0.0036543(x-2)(x-2)(x-2) x属于区间[3,4] 时；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)(x-3)-0.0409245(x-3)(x-3)(x-3) x属于区间[4,5] 时；S(x)=5.22+0.386786(x-4)-0.174138(x-4)(x-4)+0.107352(x-4)(x-4)(x-4) x属于区间[5,6] 时；S(x)=5.54+0.360567(x-5)+0.147919(x-5)(x-5)-0.268485(x-5)(x-5)(x-5) x属于区间[6,7] 时；S(x)=5.78-0.149051(x-6)-0.657537(x-6)(x-6)+0.426588(x-6)(x-6)(x-6) x属于区间[7,8] 时；S(x)=5.4-0.184361(x-7)+0.622227(x-7)(x-7)-0.267865(x-7)(x-7)(x-7) x属于区间[8,9] 时；S(x)=5.57+0.256496(x-8)-0.181369(x-8)(x-8)+0.0548728(x-8)(x-8)(x-8) x属于区间[9,10] 时；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)(x-9)+0.0583752(x-9)(x-9)(x-9) S(0.5)=2.90856 S(1.5)=3.67843 S (2.5)=4.38147 S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237 S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976 S(9.5)=5.7323

东南大学数值分析上机实验题(下)

第六章

21．（上机题）常微分方程初值问题数值解 （1） …… 此处隐藏：6422字，全部文档内容请下载后查看。喜欢就下载吧 ……