?? shuzhi2.cpp
字號:
#include "stdio.h"
#include "math.h"
void main()
{
double x0[14]={0.00,0.00,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00};/*原始數據x*/
double y[14]={0.00,0.00,5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00};/*原始數據y*/
double xx1[200]; /*定義數組x 用此存放按單位變化的x*/
double yy1[200],yy2[200],yy3[200],yy4[200]; /*定義數組y,用來存放相應的y*/
double x,px,middle,h0,px1,px2,px3,px4,hi,hi1;
int i,j,k,n;
double a[13][13]; /*用于Doolittle分解法求解方程組*/
double alfa[13],gma[13],beta[13],l[13],m0[13],p[13],q[13],r[13],s[13]; /*用于三次樣條插值法*/
/*********************************************分段線性插值**************************************************/
for (i=1;i<192;i++)
{ xx1[i]=i-1;}
k=1;
for (x=0;x<191;x++)
{
for (i=1;i<14;i++) /*確定x 屬于哪一段劃分區間*/
{ if (x<x0[i])
{n=i-1;
break;}
}
px=((x-x0[n+1])*y[n]/(x0[n]-x0[n+1]))+((x-x0[n])*y[n+1]/(x0[n+1]-x0[n])); /*求出相應區間對應的一次插值多項式 */
yy1[k]=px; /*得到y值*/
k++;
}
yy1[k]=0.00;
/*for (i=1;i<192;i++)
{/*printf("%.2f,%.2f ",xx1[i],yy1[i]);
printf("%.2f ",yy1[i]);
if (i%10==0) printf("\n");
}*/
/*********************************分段二次多項式插值*************************************/
/*for (i=2;i<14;i++)
{x1[i-1]=(x0[i-1]+x0[i])/2;
y1[i-1]=(y[i-1]+y[i])/2;}
k=1;
for (x=0;x<191;x++)
{
for (i=2;i<12;i++)
{
if (x<=x1[i])
{n=i-1;
break;}
else n=11;
}
px1=((x-x0[n+1])*(x-x1[n+1])*y1[n]/((x1[n]-x0[n+1])*(x1[n]-x1[n+1])));
px2=((x-x1[n])*(x-x1[n+1])*y[n+1]/((x0[n+1]-x1[n])*(x0[n+1]-x1[n+1])));
px3=((x-x1[n])*(x-x0[n+1])*y1[n+1]/((x1[n+1]-x1[n])*(x1[n+1]-x0[n+1])));
px=px1+px2+px3;
yy2[k]=px;
k++;
}
/* for (i=1;i<192;i++)
{printf("%.2f,%.2f ",xx1[i],yy2[i]);}
*/
k=1;
for (x=0;x<191;x++)
{
for (i=2;i<14;i++)
{
if (x<=(x0[1]+x0[2])/2) n=1;
else if ((x<=(x0[i]+x0[i+1])/2)&(x>(x0[i-1]+x0[i])/2)) /*確定x 屬于哪一段劃分區間*/
{n=i-1;
break;}
else n=11;
}
px1=(x-x0[n+1])*(x-x0[n+2])*y[n]/((x0[n]-x0[n+1])*(x0[n]-x0[n+2]));
px2=(x-x0[n])*(x-x0[n+2])*y[n+1]/((x0[n+1]-x0[n])*(x0[n+1]-x0[n+2]));
px3=(x-x0[n])*(x-x0[n+1])*y[n+2]/((x0[n+2]-x0[n])*(x0[n+2]-x0[n+1]));
px=px1+px2+px3; /*求出相應區間對應的二次插值多項式 */
yy2[k]=px; /*得到y值*/
k++;
}
/*for (i=1;i<192;i++)
{/*printf("(%.2f,%.2f) ",xx1[i],yy2[i]);
printf("%.2f ",yy2[i]);}*/
/********************************分段三次插值**********************************************/
k=1;
for (x=0;x<191;x++)
{
for (i=3;i<13;i++)
{ if (x<x0[i])
{n=i-2;
break;}
else n=10; /*確定x 屬于哪一段劃分區間*/
}
px1=(x-x0[n+1])*(x-x0[n+2])*(x-x0[n+3])*y[n]/((x0[n]-x0[n+1])*(x0[n]-x0[n+2])*(x0[n]-x0[n+3]));
px2=(x-x0[n])*(x-x0[n+2])*(x-x0[n+3])*y[n+1]/((x0[n+1]-x0[n])*(x0[n+1]-x0[n+2])*(x0[n+1]-x0[n+3]));
px3=(x-x0[n])*(x-x0[n+1])*(x-x0[n+3])*y[n+2]/((x0[n+2]-x0[n])*(x0[n+2]-x0[n+1])*(x0[n+2]-x0[n+3]));
px4=(x-x0[n])*(x-x0[n+1])*(x-x0[n+2])*y[n+3]/((x0[n+3]-x0[n])*(x0[n+3]-x0[n+1])*(x0[n+3]-x0[n+2]));
px=px1+px2+px3+px4; /*求出相應區間對應的二次插值多項式 */
yy3[k]=px; /*得到y值*/
k++;
}
/*for (i=1;i<192;i++)
{printf("%.2f ",yy3[i]);}*/
/******************************三次樣條插值************************************************/
for (i=0;i<13;i++)
{
for (j=0;j<13;j++)
{a[i][j]=0.0;} /*初始化系數矩陣*/
}
for (i=1;i<12;i++)
{hi=x0[i+1]-x0[i];
hi1=x0[i+2]-x0[i+1];
alfa[i]=hi1/(hi+hi1);
gma[i]=1-alfa[i];
beta[i]=(6/(hi+hi1))*(((y[i+2]-y[i+1])/hi1)-(y[i+1]-y[i])/hi); /*求出矩陣中的alfa,beta,gma值*/
}
alfa[12]=(x0[2]-x0[1])/(x0[2]-x0[1]+x0[13]-x0[12]);
gma[12]=1-alfa[12];
beta[12]=(6/(x0[2]-x0[1]+x0[13]-x0[12]))*((y[2]-y[1])/(x0[2]-x0[1])-(y[13]-y[12])/(x0[13]-x0[12]));
for (i=1;i<13;i++)
{
a[i][i+1]=alfa[i];
a[i][i-1]=gma[i];
for (j=1;j<13;j++)
{
if (i==j) {a[i][j]=2;} /*得到系數矩陣和右端項*/
}
}
a[1][12]=gma[1];
a[12][1]=alfa[12];
/*for (i=1;i<13;i++)
{
for (j=1;j<13;j++)
{printf("%f ",a[i][j]);}
printf("\n");
}
for (i=1;i<13;i++)
{printf("%f ",beta[i]);}
printf("\n");*/
/*********************利用Doolittle分解法反解方程求得新的u0********************************************/
/******分解過程***********/
/*for (j=1;j<13;j++)
{
u1[1][j]=a[1][j];
}
for (i=2;i<13;i++)
{
l1[i][1]=a[i][1]/u1[1][1];
}
for (im=1;im<13;im++)
{
for (j=im;j<13;j++)
{
middle=0.0;
for (tt=1;tt<im;tt++)
{ middle=middle+l1[im][tt]*u1[tt][j];}
u1[im][j]=a[im][j]-middle;
}
if (im<12)
{
for (i=im+1;i<13;i++)
{
middle=0.0;
for (tt=1;tt<im;tt++)
{
middle=middle+l1[i][tt]*u1[tt][im];}
l1[i][im]=(a[i][im]-middle)/u1[im][im];
}
}
}*/
/*************分解過程結束************************************************/
/*****求解過程*****/
/*l[1]=beta[1];
for (i=2;i<13;i++)
{
middle=0.0;
for (tt=1;tt<i;tt++)
{middle=middle+l1[i][tt]*l[tt];}
l[i]=beta[i]-middle;}
m0[12]=l[12]/u1[12][12];
for (i=11;i>0;i--)
{
middle=0.0;
for (tt=i+1;tt<13;tt++)
{middle=middle+u1[i][tt]*m0[tt];}
m0[i]=(l[i]-middle)/u1[i][i];
}
m0[0]=m0[12];
for (i=0;i<13;i++)
{printf("%f ",m0[i]);
}
printf("\n");*/
n=12;
p[1]=2;
for (i=1;i<n-1;i++)
{ q[i]=alfa[i]/p[i];
p[i+1]=2-gma[i+1]*q[i];}
s[1]=gma[1]/p[1];
for (i=2;i<n-1;i++)
{ s[i]=-gma[i]*s[i-1]/p[i];
}
s[11]=(alfa[11]-gma[11]*s[10])/p[11];
r[1]=alfa[12];
for (j=2;j<n-1;j++)
{ r[j]=-r[j-1]*q[j-1];}
r[n-1]=gma[n-1]-r[n-2]*q[n-2];
middle=0.0;
for (j=1;j<n;j++)
{ middle=middle+r[j]*s[j];}
r[n]=2-middle;
l[1]=beta[1]/p[1];
for (i=2;i<n;i++)
{ l[i]=(beta[i]-gma[i]*l[i-1])/p[i];
}
middle=0.0;
for (j=1;j<n;j++)
{ middle=middle+r[j]*l[j];}
l[n]=(beta[n]-middle)/r[n];
m0[n]=l[n];
m0[n-1]=l[n-1]-s[n-1]*m0[n];
for (i=n-2;i>0;i--)
{ m0[i]=l[i]-q[i]*m0[i+1]-s[i]*m0[n];}
m0[0]=m0[12];
for (i=0;i<13;i++)
{printf("%f ",m0[i]);
}
printf("\n");
/****求解過程結束，得到m0************/
k=1;
for (x=0;x<191;x++)
{
for (i=1;i<14;i++)
{
if (x<x0[i])
{n=i-1;
break;}
}
h0=x0[n+1]-x0[n];
/*px=m0[n]*(x0[n+1]-x)*(x0[n+1]-x)*(x0[n+1]-x)/(6*h0)+m0[n+1]*(x-x0[n])*(x-x0[n])*(x-x0[n])/(6*h0)+(y[n]/h0-m0[n]*h0/6)*(x0[n+1]-x)+(y[n+1]/h0-m0[n+1]*h0/6)*(x-x0[n]);/*得到s*/
px=m0[n-1]*(x0[n+1]-x)*(x0[n+1]-x)*(x0[n+1]-x)/(6*h0)+m0[n]*(x-x0[n])*(x-x0[n])*(x-x0[n])/(6*h0)+(y[n]/h0-m0[n-1]*h0/6)*(x0[n+1]-x)+(y[n+1]/h0-m0[n]*h0/6)*(x-x0[n]);
yy4[k]=px;/*得到y*/
k++;
}
for (i=1;i<192;i++)
{printf("%.2f ",yy4[i]);
if (i%10==0) printf("\n");}
}
/*******************************主程序結束*******************/?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -