給你一段代碼,已驗(yàn)證通過!
成都創(chuàng)新互聯(lián)公司成立10余年來,這條路我們正越走越好,積累了技術(shù)與客戶資源,形成了良好的口碑。為客戶提供網(wǎng)站設(shè)計(jì)、成都網(wǎng)站制作、網(wǎng)站策劃、網(wǎng)頁設(shè)計(jì)、域名注冊、網(wǎng)絡(luò)營銷、VI設(shè)計(jì)、網(wǎng)站改版、漏洞修補(bǔ)等服務(wù)。網(wǎng)站是否美觀、功能強(qiáng)大、用戶體驗(yàn)好、性價比高、打開快等等,這些對于網(wǎng)站建設(shè)都非常重要,成都創(chuàng)新互聯(lián)公司通過對建站技術(shù)性的掌握、對創(chuàng)意設(shè)計(jì)的研究為客戶提供一站式互聯(lián)網(wǎng)解決方案,攜手廣大客戶,共同發(fā)展進(jìn)步。
你可以按照其中的算法進(jìn)行相應(yīng)的修改。
x=[0.40,0.55,0.65,0.80,0.90,1.05];%插值節(jié)點(diǎn)
y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382];%插值
n=length(x);%節(jié)點(diǎn)個數(shù)
z=input('請輸入需要計(jì)算的函數(shù)值點(diǎn)z=');
newton=[x',y'];%給出差商表的前兩列
for j=2:1:n
for i=n:-1:j
y(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));
end
end
u=y(n);
for i=n-1:-1:1
u=y(i)+u*(z-x(i));
end
fprintf('f(%8.5f)=%8.5f\n',z,u)
double newton(double *x, double *y, int n, double num, int cur, int pointNum, double answer)
{
//計(jì)算均差
for(int i = pointNum -1; icur; i--)
{
y[i] = ( y[i] - y[i-1] ) / ( x[i] - x[i-1] );
}
//已經(jīng)計(jì)算完cur自加
cur++;
//temp進(jìn)行臨時計(jì)算
double temp = y[cur];
for(int i= 0; icur; i++)
{
temp *= ( num - x[i] );
}
//將臨時計(jì)算的結(jié)果加到answer
answer += temp;
//如果得到想要的結(jié)果就返回答案 否則繼續(xù)計(jì)算
if(cur==n)
return answer;
else
return newton(x,y,n,num,cur,pointNum,answer); /// 少了 return
}
自編的,都弄上來了,缺樣條插值。這里僅是函數(shù),什么控件的編程你自己弄,那實(shí)在太簡單了。
Dim aa As Double, bb As Double '分別接收findway有根區(qū)間兩端值的變量
Dim x(1) As Double '分別接收ercigenway的根
'1.0 ercigenway 求二次方程實(shí)根 -已測試
Private Sub ercigenway(a As Single, b As Single, c As Single) 'a、b、c對應(yīng)為二次方程的系數(shù)
Dim d As Double
d = b ^ 2 - 4 * a * c
If d 0 Then
MsgBox "Δ小于0,沒有實(shí)根", , "消息"
x(0) = 0: x(1) = 0
ElseIf d = 0 Then
x(0) = -b / (2 * a): x(1) = x(0)
Else
x(0) = (-b - Sgn(b) * Sqr(d)) / (2 * a): x(1) = c / (a * x(0))
End If
End Sub
'2.1 findway 等步長掃描有根區(qū)間 -已測試
Private Sub findway(ByVal a As Single, ByVal b As Single, h As Double) 'a、b分別為待掃描區(qū)間端點(diǎn),h為步長
Dim a1 As Double
a1 = a
Do
If f(a1) * f(a1 + h) = 0 Then
aa = a1: bb = a1 + h
Exit Sub
End If
a1 = a1 + h
Loop While a1 b
If a1 b Then
MsgBox "沒有找到有根區(qū)間,請換更小的步長試一下"
Exit Sub
End If
End Sub
'2.2 erfenfun 二分法求根 -已測試
Private Function erfenfun(ByVal a As Single, ByVal b As Single, eps As Double) 'a、b為有根區(qū)間端點(diǎn),eps為誤差
Dim x0 As Double, x1 As Double, x2 As Double, f0 As Double, f1 As Double, f2 As Double
x1 = a: x2 = b
Do
x0 = (x1 + x2) / 2
f0 = f(x0)
If f0 = 0 Then
Exit Do
Else
f1 = f(x1): f2 = f(x2)
If f0 * f1 0 Then
x2 = x0
Else
x1 = x0
End If
End If
Loop While Abs(x1 - x2) eps
x0 = (x1 + x2) / 2
erfenfun = x0
End Function
'2.4 newtonfxfun Newton切線法 -已測試
Private Function newtonfxfun(ByVal x0 As Double, eps As Double) As Double 'x0為附近根,eps為誤差
Dim x1 As Double, f0 As Double, f1 As Double
x1 = x0
Do
x0 = x1
f0 = f(x0): f1 = fd(x0) 'fd表示f的導(dǎo)函數(shù)
If Abs(f1) eps Then
x1 = x0: Exit Do
End If
x1 = x0 - f0 / f1
Loop Until Abs(x1 - x0) eps
newtonfxfun = x1
End Function
'2.3 stediedaifun Seffensen加速迭代法 (方程形式為x-f(x)=0) -已測試
Private Function stediedaifun(ByVal x0 As Double, eps1 As Double, eps2 As Double) As Double 'x0為解析解附近的根,eps1為輸出結(jié)果誤差,eps2為迭代能否繼續(xù)判斷標(biāo)準(zhǔn)
Dim y As Double, z As Double, x1 As Double
x1 = x0
Do
x0 = x1
y = f(x0): z = f(y)
If Abs(z - 2 * y + x0) eps2 Then
MsgBox "為滿足eps2條件,不能繼續(xù)迭代"
Exit Function
End If
x1 = x0 - (y - x0) ^ 2 / (z - 2 * y + x0)
Loop Until Abs(x1 - x0) eps1
stediedaifun = x1
End Function
'2.5 newtonfxnfun n次代數(shù)方程N(yùn)ewton切線法 -已測試
Private Function newtonfxnfun(a() As Single, eps As Double, x0 As Double) As Double 'a()分別存儲按降冪排列的方程的n個系數(shù),eps為誤差,x0為附近根
Dim k As Integer, n As Integer, f0 As Double, f1 As Double, x1 As Double
n = UBound(a)
x1 = x0
Do
x0 = x1
f0 = a(0): f1 = f0
For k = 1 To n - 1
f0 = a(k) + f0 * x0
f1 = f0 + f1 * x0
Next k
f0 = a(n) + f0 * x0
x1 = x0 - f0 / f1
Loop Until Abs(x1 - x0) eps
newtonfxnfun = x1
End Function
'2.6 linecutfun 弦截法 -已測試
Private Function linecutfun(ByVal x0 As Double, ByVal x1 As Double, eps As Double, n As Long) As Double 'n為迭代次數(shù)限制,x0、x1為有根區(qū)間端點(diǎn),eps為誤差
Dim f0 As Double, f1 As Double, f2 As Double
Dim x2 As Double, i As Long
f0 = f(x0): f1 = f(x1)
For i = 1 To n
x2 = x1 - (x1 - x0) * f1 / (f1 - f0)
f2 = f(x2)
If Abs(f2) eps Then
Exit For
End If
x0 = x1: x1 = x2: f0 = f1: f1 = f2
Next i
If i = n + 1 Then
MsgBox "要求的計(jì)算次數(shù)太低,沒有達(dá)到精度要求"
End If
linecutfun = x2
End Function
'4.1 lagrangeczfun 拉格朗日插值法 -已測試
Private Function lagrangeczfun(a() As Double, ByVal u As Double) As Double 'a(1,n)存儲n+1個節(jié)點(diǎn),u為插值點(diǎn)
Dim i As Integer, j As Integer, n As Integer
Dim l As Double, v As Double
v = 0
n = UBound(a, 2)
For j = 0 To n
l = 1#
For i = 0 To n
If i = j Then GoTo hulue
l = l * (u - a(0, i)) / (a(0, j) - a(0, i))
hulue:
Next i
v = v + l * a(1, j)
Next j
lagrangeczfun = v
End Function
'4.2 newtonczfun newton插值法 -已測試
Private Function newtonczfun(a() As Double, u As Double) As Double 'a(1,n)存儲n+1個節(jié)點(diǎn),u為插值點(diǎn)
Dim n As Integer, i As Integer, j As Integer, k As Integer
Dim z() As Double, f() As Double, v As Double
n = UBound(a, 2)
ReDim z(n), f(n)
For i = 0 To n
z(i) = a(1, i)
Next i
For i = 1 To n
k = k + 1
For j = i To n
f(j) = (z(j) - z(j - 1)) / (a(0, j) - a(0, j - k))
Next j
For j = i To n
z(j) = f(j)
Next j
Next i
f(0) = a(1, 0)
v = 0
For i = n To 0 Step -1
v = v * (u - a(0, i)) + f(i)
Next i
newtonczfun = v
End Function
'4.3 hermiteczfun Hermite插值法 -已測試
Private Function hermiteczfun(a() As Double, fd() As Double, u As Double) As Double 'a(1,n)存儲n+1個節(jié)點(diǎn),fd(n)存儲n+1個節(jié)點(diǎn)處導(dǎo)數(shù)值,u為插值點(diǎn)
Dim l() As Double, ld() As Double, g() As Double, h() As Double, aim As Double
Dim n As Integer, i As Integer, j As Integer
n = UBound(a)
ReDim l(n), ld(n), g(n), h(n)
aim = 0
For i = 0 To n
l(i) = 1: ld(i) = 0
For j = 0 To n
If j = i Then GoTo hulue
l(i) = l(i) * (u - a(0, j)) / (a(0, i) - a(0, j))
ld(i) = ld(i) + 1 / (a(0, i) - a(0, j))
hulue:
Next j
g(i) = (1 + 2 * (a(0, i) - u) * ld(i)) * l(i) * l(i)
h(i) = (u - a(0, i)) * l(i) * l(i)
aim = aim + g(i) * a(1, i) + h(i) * fd(i)
Next i
hermiteczfun = aim
End Function
'5.2.1 tixingjffun 變步長梯形積分法 -已測試
Private Function tixingjffun(a As Single, b As Single, eps As Double, m As Long) As Double 'a、b分別為積分上下限,eps為誤差,m為最大計(jì)算次數(shù)
Dim h As Double, t1 As Double, t2 As Double, t As Double, hh As Double
Dim n As Long: n = 1
h = b - a: t1 = h * (f(a) + f(b)) / 2
Do
t = 0
For i = 1 To n
t = t + f(a + (i - 0.5) * h)
Next i
hh = h * t
t2 = (t1 + hh) / 2
If Abs(t2 - t1) eps Then Exit Do
t1 = t2: h = h / 2: n = 2 * n
Loop Until n 2 * m
If n 2 * m Then
MsgBox "計(jì)算次數(shù)預(yù)定太小,不能達(dá)到誤差要求"
End If
tixingjffun = t2
End Function
'5.2.2 simpsonjffun 變步長Simpson積分法 -已測試
Private Function simpsonjffun(a As Single, b As Single, eps As Double, m As Long) As Double 'a、b分別為積分上下限,eps為誤差,m為最大計(jì)算次數(shù)
Dim n As Long, i As Long
Dim h As Double, t1 As Double, t2 As Double, hh As Double, s1 As Double, s2 As Double
n = 1: h = b - a: t1 = h * (f(a) + f(b)) / 2
hh = h * (f((a + b) / 2)): s1 = (t1 + 2 * hh) / 3
Do
n = 2 * n: h = h / 2: t2 = (t1 + hh) / 2
t = 0
For i = 1 To n
t = t + f(a + (i - 0.5) * h)
Next i
hh = t * h
s2 = (t1 + 2 * hh) / 3
If Abs(s2 - s1) eps Then Exit Do
t1 = t2: s1 = s2
Loop Until n m
If n m Then MsgBox "計(jì)算次數(shù)預(yù)定太小,不能達(dá)到誤差要求"
simpsonjffun = s2
End Function
'5.3 Rombergjffun Romberg積分法
Private Function rombergjffun(a As Single, b As Single, eps As Double) As Double
Dim k As Integer, n As Integer, h As Double
k = 0: n = 1: h = b - a
End Function
'5.5.1 ds1fun 求一階導(dǎo)數(shù) -已測試
Private Function ds1fun(x0 As Single, eps As Double) As Double 'x0為求導(dǎo)點(diǎn),eps為誤差
Dim h As Double, t1 As Double, t2 As Double
h = 1: t1 = (f(x0 + h) - f(x0 - h)) / (2 * h)
h = h / 2: t2 = (f(x0 + h) - f(x0 - h)) / (2 * h)
Do While Abs(t2 - t1) eps
t1 = t2
h = h / 2
t2 = (f(x0 + h) - f(x0 - h)) / (2 * h)
Loop
ds1fun = t2
End Function
'5.5.2 ds2fun 求二階導(dǎo)數(shù) -已測試
Private Function ds2fun(x0 As Single, eps As Double) As Double 'x0為求導(dǎo)點(diǎn),eps為誤差
Dim h As Double, t1 As Double, t2 As Double
h = 1: t1 = (f(x0 + h) + f(x0 - h) - 2 * f(x0)) / (h * h)
h = h / 2: t2 = (f(x0 + h) + f(x0 - h) - 2 * f(x0)) / (h * h)
Do While Abs(t2 - t1) eps
t1 = t2
h = h / 2
t2 = (f(x0 + h) + f(x0 - h) - 2 * f(x0)) / (h * h)
Loop
ds2fun = t2
End Function
牛頓插值法:
#includestdio.h
#includealloc.h
float Language(float *x,float *y,float xx,int n)
{
int i,j;
float *a,yy=0.0;
a=(float *)malloc(n*sizeof(float));
for(i=0;i=n-1;i++)
{
a[i]=y[i];
for(j=0;j=n-1;j++)
if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);
yy+=a[i];
}
free(a);
return yy;
}
void main()
{
float x[4]={0.56160,0.5628,0.56401,0.56521};
float y[4]={0.82741,0.82659,0.82577,0.82495};
float xx=0.5635,yy;
float Language(float *,float *,float,int);
yy=Language(x,y,xx,4);
printf("x=%f,y=%f\n",xx,yy);
getchar();
}
2.牛頓插值法#includestdio.h
#includemath.h
#define N 4
void Difference(float *x,float *y,int n)
{
float *f;
int k,i;
f=(float *)malloc(n*sizeof(float));
for(k=1;k=n;k++)
{
f[0]=y[k];
for(i=0;ik;i++)
f[i+1]=(f[i]-y[i])/(x[k]-x[i]);
y[k]=f[k];
}
return;
}
main()
{
int i;
float varx=0.895,b;
float x[N+1]={0.4,0.55,0.65,0.8,0.9};
float y[N+1]={0.41075,0.57815,0.69675,0.88811,1.02652};
Difference(x,(float *)y,N);
b=y[N];
for(i=N-1;i=0;i--)b=b*(varx-x[i])+y[i];
printf("Nn(%f)=%f",varx,b);
getchar();
}
留下個郵箱,我發(fā)給你:牛頓插值法的程序設(shè)計(jì)與應(yīng)用
class?Lagrange
{
public?static?double?esitimate(double?x,int?n,double?k[][])
{
double?y=0;
int?i;
int?j;
double?w;
double?w1;
for(i=0;in;i++)
{
w=1;
w1=1;
for(j=0;jn;j++)
{
if(j==i)
continue;
else
{
w=w*(x-k[j][0]);
w1=w1*(k[i][0]-k[j][0]);
}
}
//System.out.printf("y=%f,w=%f,w1=%f\n",k[i][1],w,w1);
y=y+k[i][1]*w/w1;
}
return?y;
}
}
class?Newton
{
double?ad[][];//均差表
int?n;
Newton(int?n,double?k[][])
{
int?i,j;
this.n=n;
ad=new?double[n][n];
for(i=0;in;i++)
ad[i][0]=k[i][1];
for(i=1;in;i++)
ad[i][1]=(ad[i][0]-ad[i-1][0])/(k[i][0]-k[i-1][0]);
for(i=2;in;i++)
for(j=i;jn;j++)
{
ad[j][i]=(ad[j][i-1]-ad[j-1][i-1])/(k[j][0]-k[j-i][0]);
//System.out.printf("%f?%f?%f?%f\n",?ad[j][i-1],ad[j-1][i-1],k[j][0],k[j-i][0]);
}
show();
}
void?update(double?x,double?y,double?k[][])
{
int?i,j;
n++;
k[n-1][0]=x;
k[n-1][1]=y;
double?temp[][]=new?double[n][n];
for(i=0;in-1;i++)
for(j=0;j=i;j++)
temp[i][j]=ad[i][j];
temp[n-1][0]=y;
for(i=1;in;i++)
temp[n-1][i]=(temp[n-1][i-1]-temp[n-2][i-1])/(x-k[n-1-i][0]);
ad=temp;
show();
}
void?show()
{
int?i,j;
System.out.println("均差計(jì)算可列均差表如下:");
for(i=0;in;i++)
{
for(j=0;j=i;j++)
System.out.printf("%f???",?ad[i][j]);
System.out.printf("\n");
}
System.out.println("*******************************************");
}
double?esitimate(double?x,double?k[][])
{
int?i,j;
double?y=ad[0][0];
double?w;
for(i=1;in;i++)
{
w=ad[i][i];
for(j=0;ji;j++)
w=w*(x-k[j][0]);
y=y+w;
}
return?y;??
}
}
public?class?interpolate?
{
public?static?void?main(String?args[])
{
int?n=20;
double?k[][]=new?double[n][2];
int?i;
n=3;
k[0][0]=0.32;
k[0][1]=0.314567;
k[1][0]=0.34;
k[1][1]=0.333487;
k[2][0]=0.36;
k[2][1]=0.352274;
System.out.println("拉格朗日插值的節(jié)點(diǎn):");
for(i=0;in;i++)
System.out.printf("%f????%f\n",?k[i][0],k[i][1]);
System.out.println("估算0.3367處的函數(shù)值:");
n=2;
System.out.printf("采用%d次插值得:%f\n",n-1,Lagrange.esitimate(0.3367,n,k));
n=3;
System.out.printf("采用%d次插值得:%f\n",n-1,Lagrange.esitimate(0.3367,n,k));
System.out.println("*******************************************");
k[0][0]=0.4;
k[0][1]=0.41075;
k[1][0]=0.55;
k[1][1]=0.57815;
k[2][0]=0.65;
k[2][1]=0.69675;
k[3][0]=0.80;
k[3][1]=0.88811;
k[4][0]=0.90;
k[4][1]=1.02652;
k[5][0]=1.05;
k[5][1]=1.25382;
Newton?nt=new?Newton(3,k);
nt.update(0.80,?0.88811,?k);
nt.update(0.90,?1.02652,?k);
nt.update(1.05,?1.25382,?k);
nt=new?Newton(5,k);
System.out.println("估算0.596處的函數(shù)值:");
System.out.printf("%f",nt.esitimate(0.596,?k));
}
}