第3章 函数近似方法(附录)

第3章附录
3.1.5 n+1点n次插值
function y = lagrange_2(x0,y0,x)
n = length(x0); m = length(x);
for i = 1:m
z = x(i); s = 0.0;
for k = 1:n
p = 1.0;
for j =1:n
if j~=k
p = p*(z-x0(j))/(x0(k)-x0(j));
end
end
s = s+ p*y0(k);
end
y(i) = s;
End
============================================== 3.1.6 三次样条插值
例题3.1.7计算程序
% demo_spline.m
x = -1:0.1:1;
y = 1./(1+x.*x);
xx = -1:.01:1;
yy = spline(x,y,xx);
plot(x,y,'bo',xx,yy,'r','LineWidth',2)
title('y=1/(1+x^2): 样条插值…
','FontSize',12);
xlabel('x','FontSize',12);
ylabel('y','FontSize',12)
% demo_splines.m
x = [-1.:0.1:1.];
n = length(x);
y = 1./(1.+x.*x);
m = 10*(n-1)+1;
dx = (x(n)-x(1))/(m-1.)
for i = 2:n-1
a(i) = x(i+1)-x(i);
b(i) = 2.*(x(i+1)-x(i));
c(i) = x(i+1)-x(i);
f(i) = 6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i+1)-x(i)));
end
a(1) = 0.; b(1) = 1.;c(1) = 0.; f(1) = 0.;
a(n) = 0.; b(n) = 1.;c(n) = 0.; f(n) = 0.;
ypp = tri(a,b,c,f);
for i = 1:n-1
c1(i) = (ypp(i+1)-ypp(i))/(x(i+1)-x(i));
c2(i) = (x(i+1)*ypp(i)-x(i)*ypp(i+1))/(x(i+1)-x(i));
c3(i) = ((x(i)*x(i)-2.*x(i+1)*x(i+1)-2.*x(i)*x(i+1))*ypp(i)+(2.*x(i)*x(i)-x(i+1)*x(i+1)...
+ 2.*x(i)*x(i+1))*ypp(i+1)+6.*(y(i+1)-y(i)))/6./(x(i+1)-x(i));
c4(i) = -(x(i)*x(i+1)*((x(i)-2.*x(i+1))*ypp(i)+(2.*x(i)-x(i+1))*ypp(i+1))...
+ 6.*(x(i)*y(i+1)-x(i+1)*y(i)))/6./(x(i+1)-x(i));
end
j =1;
for i = 1:n-1
dj = floor((x(i+1)-x(i)+dx/2)/dx);
for k =j:(j+dj-1)
z(k) = x(1)+(k-1)*dx;
s(k) = c1(i)*z(k)*z(k)*z(k)/6.+c2(i)*z(k)*z(k)/2.+c3(i)*z(k)+c4(i);
sp(k)= c1(i)*z(k)*z(k)/2.+c2(i)*z(k)+c3(i);
spp(k)=c1(i)*z(k)+c2(i);
end
j = j+dj;
end
yy = 1./(1+z.*z);
plot(z,yy,'k',x,y,'bo',z,s,'r','LineWidth',2)
xlabel('x','FontSize',12);
ylabel('y','FontSize',12)
legend('原曲线','插值点','插值曲线');
title('y=1/(1+x^2):样条插值','FontSize',12);
==============================================
例题3.2.1计算程序
!!!!linear_fit.for!!!
program linear_fit
! linear_fit
implicit none
integer n,i,m
parameter (n=7,m=20)
real, dimension (n):: x,y
real, dimension (m):: xi,yi
real a,b,dx
open(5,file='exa3_2_1.dat')
open(2,file='exa3_2_1_old.dat')
data x/0.5, 1.2, 2.1, 2.9, 3.6, 4.5, 5.7/ data y/2.81,3.24,3.80,4.30,4.73,5.29,6.03/ write(5,'(2f14.6)') (x(i),y(i),i=1,n)
call fit(n,x,y,a,b)
print *,'a=',a,' b=',b
dx=(x(n)-x(1))/(m-1)
do i=1,m
xi(i)=x(1)+(i-1)*dx
yi(i)=a+b*xi(i)
end do
write(2,'(2f14.6)') (xi(i),yi(i),i=1,m) end
subroutine fit(n,x,y,a,b)
implicit none
integer n,i
real, dimension (n):: x,y
real xp,yp,x2,xyp,a,b,c
xp = 0.0
yp = 0.0
x2 = 0.0
xyp= 0.0
do i = 1,n
x2 = x2+x(i)*x(i)
xp = xp+x(i)
yp = yp+y(i)
xyp= xyp+x(i)*y(i)
end do
xp = xp/n
yp = yp/n
x2 = x2/n
xyp= xyp/n
c = x2-xp*xp
a = (yp*x2-xp*xyp)/c
b = (xyp-xp*yp)/c
return
end
//*linear fit.cpp*/
#include <stdio.h>
main()
{
int n =7,i;
float x[]={0.5, 1.2, 2.1, 2.9, 3.6, 4.5, 5.7};
float y[]={2.81,3.24,3.80,4.30,4.73,5.29,6.03};
float sx=0.,sy=0.,sxy=0.,sx2=0,deno,a,b;
for(i=0;i<=6;i++)
{
sx=sx+x[i];
sy=sy+y[i];
sxy=sxy+x[i]*y[i];
sx2=sx2+x[i]*x[i];
}
deno=n*sx2-sx*sx;
a=(sy*sx2-sx*sxy)/deno;
b=(n*sxy-sy*sx)/deno;
printf("a=%6.2f b=%6.2f\n",a,b);
return(0);
}
============================================== 例题3.2.2计算程序
! nonlinearfit.for
program nonlinearfit
! non-linear_fit
implicit none
integer n,i,m
parameter (n=8,m=20)
real, dimension (n):: x,y,yc
real a,b,dx
data x/1,2,3,4,5,6,7,8/
data y/15.3,20.5,27.4,36.6,49.1,65.6,87.8,117.6/
do i=1,n
yc(i)=log(y(i))
end do
call fit(n,x,yc,a,b)
a=exp(a)
print *,'a=',a,' b=',b
end
subroutine fit(n,x,y,a,b)
implicit none
integer n,i
real, dimension (n):: x,y
real xp,yp,x2,xyp,a,b,c
xp = 0.0
yp = 0.0
x2 = 0.0
xyp= 0.0
do i = 1,n
x2 = x2+x(i)*x(i)
xp = xp+x(i)
yp = yp+y(i)
xyp= xyp+x(i)*y(i)
end do
xp = xp/n
yp = yp/n
x2 = x2/n
xyp= xyp/n
c = x2-xp*xp
a = (yp*x2-xp*xyp)/c
b = (xyp-xp*yp)/c
return
end
============================================= 3.2.3 m次多项式拟合
例题3.2.3计算程序
!!! polyfit.For
program polyfit
! 多项式拟合
implicit none
integer n,m,i,l,k,j
parameter (n=11,m=6,k=21)
real, dimension (n):: x,y
real, dimension (k):: xn,yn,xo,yo
real s(m,m),t(m),js(m),z(m),dx
open(5,file='exa3_2_3.dat')
open(2,file='exa3_2_3_old.dat')
data x/-1.0,-0.8,-0.6,-0.4,-0.2,0.0,0.2,0.4,0.6,0.8,1.0/ do i=1,n
y(i)=exp(x(i))
end do
call fitp(n,m,x,y,s,t)
call gaus(s,t,m,z,l,js)
if (l.ne.0) then
write(*,10) (i,z(i),i=1,6)
end if
10 format(1x,'x(',i2,' )=',d15.6)
dx=2.0/(k-1)
do i=1,k
xn(i)=-1.0+dx*(i-1)
xo(i)=xn(i)
yo(i)=exp(xo(i))
yn(i)=z(m)
do j=1,m-1
yn(i)=yn(i)*xn(i)+z(m-j)
end do
write(5,'(2f14.6)') xn(i),yn(i)
write(2,'(2f14.6)') xo(i),yo(i) end do
end
subroutine fitp(n,m1,x,y,s,t)
implicit none
integer n,m1,m,j,i,i1,mi
real, dimension (n):: x,y
real s(m1,m1),t(m1)
m=m1-1
s=0.0
t=0.0
s(1,1)=n
do j=1,n
t(1)=t(1)+y(j)
end do
do i=2,m1
i1=i-1
mi=m1+i-2
do j=1,n
s(i,1)=s(i,1)+x(j)**i1
s(m1,i)=s(m1,i)+x(j)**mi
t(i)=t(i)+x(j)**i1*y(j)
end do
end do
do j=2,m
do i=j,m
s(i,j)=s(i+1,j-1)
end do
end do
do i=1,m
i1=i+1
do j=i1,m1
s(i,j)=s(j,i)
end do
end do
return
end
subroutine gaus(a,b,n,x,l,js)
real a(n,n),x(n),b(n),js(n)
real t
l=1
do 50 k=1,n-1
d=0.0
do 210 i=k,n
do 210 j=k,n
if (abs(a(i,j)).gt.d) then
d=abs(a(i,j))
js(k)=j
is=i
end if
210 continue
if (d+1.0.eq.1.0) then
l=0
else
if (js(k).ne.k) then
do 220 i=1,n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
220 continue
end if
if (is.ne.k) then
do 230 j=k,n
t=a(k,j)
a(k,j)=a(is,j)
a(is,j)=t
230 continue
t=b(k)
b(k)=b(is)
b(is)=t
end if
end if
if (l.eq.0) then
write(*,100)
return
end if
do 10 j=k+1,n
a(k,j)=a(k,j)/a(k,k)
10 continue
b(k)=b(k)/a(k,k)
do 30 i=k+1,n
do 20 j=k+1,n
a(i,j)=a(i,j)-a(i,k)*a(k,j) 20 continue
b(i)=b(i)-a(i,k)*b(k)
30 continue
50 continue
if (abs(a(n,n))+1.0.eq.1.0) then
l=0
write(*,100)
return
end if
x(n)=b(n)/a(n,n)
do 70 i=n-1,1,-1
t=0.0
do 60 j=i+1,n
t=t+a(i,j)*x(j)
60 continue
x(i)=b(i)-t
70 continue
100 format(1x,' fail ')
js(n)=n
do 150 k=n,1,-1
if (js(k).ne.k) then
t=x(k)
x(k)=x(js(k))
x(js(k))=t
end if
150 continue
return
end
% demo_polyfits
clc;clear all;
format long;
x =-1:0.2:1;
y = exp(x);
m = 5;
a = polyfits(x,y,m)
c =[a(6) a(5) a(4) a(3) a(2) a(1)];
x1=-1:0.02:1;
y1=polyval(c,x1);
plot(x1,y1,'r',x,y,'o');
title('5次多项式拟合exp(x)','FontSize',14); xlabel('x','FontSize',14);
ylabel('y','FontSize',14);
function a =polyfits(x,y,m)
n=length(x); d(1:(2*m+1))=0;
b(1:(m+1))=0;
for j=1:(2*m+1)
for k=1:n
d(j)=d(j)+x(k)^(j-1);
if(j<(m+2))
b(j)=b(j)+y(k)*x(k)^(j-1);
end
end
end
c(1,:)=d(1:(m+1));
for s=2:(m+1)
c(s,:)=d(s:(m+s));
end
a = gauss(c,b');。