计算方法4-6章存档
作者：互联网
第四章 插值与拟合
给定一系列的点，求多项式函数
可以用n个点确定n个未知量，但此方程组病态，误差极大。
1.插值余项（PPT1
用于估算误差（P76
2.拉格朗日插值多项式
为每个点构造一个格式统一的函数，使得取到该点时值为1，其余点均为0.
一种构造：
取2点，为线性插值，代码：
function y=xianxing(X,Y,x) x0=X(1); y0=Y(1); x1=X(2); y1=Y(2); l0=(x-x1)/(x0-x1); l1=(x-x0)/(x1-x0); y=y0*l0+y1*l1;View Code
取3点，为抛物插值，代码：
function y=paowu(X,Y,x) x0=X(1); x1=X(2); x2=X(3); y0=Y(1); y1=Y(2); y2=Y(3); L0=(x-x1)*(x-x2)/(x0-x1)/(x0-x2); L1=(x-x0)*(x-x2)/(x1-x0)/(x1-x2); L2=(x-x0)*(x-x1)/(x2-x0)/(x2-x1); y=y0*L0+y1*L1+y2*L2;View Code
误差分析：套余项式，注意代最近的点。（P79
3.分段低次插值
提高插值多项式的次数会产生龙格现象，因此换用多个低次拼凑。
分段线性插值
代码：
function yy=fdxx(x,y,xx)
n=size(x,2);
for i=1:n-1
if x(i)<xx&&xx<x(i+1)
L1=(xx-x(i+1))/(x(i)-x(i+1));
L2=(xx-x(i))/(x(i+1)-x(i));
yy=L1*y(i)+L2*y(i+1);
break;
elseif x(i)==xx
yy=y(i);
end
end
end
View Code
分段抛物插值
代码：
function y=fdpw(xi,yi,x)
n=size(xi,2);
if x<xi(2)
L1=(x-xi(1))*(x-xi(3))/(xi(1)-xi(2))/(xi(1)-xi(3));
L2=(x-xi(1))*(x-xi(3))/(xi(2)-xi(1))/(xi(2)-xi(3));
L3=(x-xi(1))*(x-xi(2))/(xi(3)-xi(1))/(xi(3)-xi(2));
y=L1*yi(1)+L2*yi(2)+L3*yi(3);
elseif x>xi(end-1)
L1=(x-xi(end-1))*(x-xi(end))/(xi(end-2)-xi(end-1))/(xi(end-2)-xi(end));
L2=(x-xi(end-2))*(x-xi(end))/(xi(end-1)-xi(end-2))/(xi(end-1)-xi(end));
L3=(x-xi(end-2))*(x-xi(end-1))/(xi(end)-xi(end-2))/(xi(end)-xi(end-1));
y=L1*yi(1)+L2*yi(2)+L3*yi(3);
else
for k=2:n-1
if xi(k+1)>x
if abs(x-xi(k))<abs(x-xi(k+1))
i=k-1;
else
i=k;
end
L1=(x-xi(i+1))*(x-xi(i+2))/(xi(i)-xi(i+1))/(xi(i)-xi(i+2));
L2=(x-xi(i))*(x-xi(i+2))/(xi(i+1)-xi(i))/(xi(i+1)-xi(i+2));
L3=(x-xi(i))*(x-xi(i+1))/(xi(i+2)-xi(i))/(xi(i+2)-xi(i+1));
y=L1*yi(i)+L2*yi(i+1)+L3*yi(i+2);
end
end
end
end
View Code
三次样条插值
分段线性和分段抛物插值会导致曲线不平滑。三次样条插值用多个三次函数拼接，保证曲线二阶导连续。一个三次函数有4个未知数，n个三次函数就有4n个未知数，当前点函数值、当前点连续、一阶导连续、二阶导连续有4n-2个已知量，还差两个，因此需要2个初值，一般为两端点的一阶导或二阶导。实际应用中用二阶导表示函数，这样只需计算n+1个未知数，从而推导出简便的式子（4.52）。（P91-94
代码：
function s=threesimple1(X,Y,x,y0,yn)
% 三次样条插值函数第一类型代码
% s函数表示三次样条插值函数插值点对应的函数值
% 根据三次样条参数函数求出的D,h,A,g,M
% x表示求解插值点函数点，X为已知插值点
[D,h,A,g,M]=Mspline3(X,Y,y0,yn)
n=length(X); m=length(x);
for t=1:m
for i=1:n-1
if (x(t)<=X(i+1))&&(x(t)>=X(i))
p1=M(i,1)*(X(i+1)-x(t))^3/(6*h(i));
p2=M(i+1,1)*(x(t)-X(i))^3/(6*h(i));
p3=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
p4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
s(t)=p1+p2+p3+p4;
break;
else
s(t)=0;
end
end
end
end
function [D,h,A,g,M]=Mspline3(X,Y,y0,yn)
% 自然边界条件的三次样条函数(第一种边界条件)
% 此函数为M值求值函数
% D,h,A,g,M输出量分别为系数矩阵D，插值宽度h，差商表A，g值,M值
n=length(X);
A=zeros(n,n);A(:,1)=Y';D=zeros(n,n);g=zeros(n,1);
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
end
end
for i=1:n-1
h(i)=X(i+1)-X(i);
end
for i=1:n
D(i,i)=2;
D(1,2)=1;
D(n,n-1)=1;
if (i==1)
g(i,1)=6/h(i)*(A(2,2)-y0);
elseif (i==n)
g(i,1)=6/h(i-1)*(yn-A(i,2));
else
g(i,1)=(6/(h(i-1)+h(i)))*(A(i+1,2)-A(i,2));
end
end
for i=1:n-2
u(i)=h(i)/(h(i)+h(i+1));
n(i)=1-u(i);
D(i+1,i+2)=n(i);
D(i+1,i)=u(i); %改到这里
end
M=D\g;
%M=[0;M;0];
end
View Code
4.曲线拟合的最小二乘法（PPT2
由于采样点有误差，因此希望得到一个更接近真实规律的近似表达式，使得偏差平方和最小。求法为先根据观察选择一个函数，再求误差表达式取得极小值的条件，即一阶导等于0推出的方程组（P104
方程组是线性的比较好解，因此我们通过换元把函数化为线性形式，例如指数形式表达式取对数等。
代码：
function p=leastsquare(x,y,m)
A=zeros(m+1,m+1);
b=zeros(1,m+1);
for i=0:m
for j=0:m
A(i+1,j+1)=sum(x.^(i+j));
end
b(i+1)=sum(x.^i.*y);
end
p=vpgauss(A,b);
View Code
第五章 数值微分与数值积分
1.数值微分（PPT1
根据离散点函数值推断某点导数近似值。
求法：通过构造插值多项式得到数值微分公式及其余项，常用其二级结论两点公式、三点公式计算导数（P117
2.数值积分（PPT2
根据离散点函数值推断区域上函数积分值。
求法：同样通过构造插值多项式代替被积函数，得到插值型求积公式及其余项（P121-122
数值积分公式的代数精度
以代数多项式为准，对任意不高于m次代数多项式准确成立，而m+1次不成立，称该数值积分公式代数精度为m。例：P123
牛顿-科茨方法
用n次插值函数多项式近似被积函数。
复合就是将积分区域分成n段，每段按对应方法积分。误差是步长的n次方，称为n阶方法。
（复合）梯形公式
用梯形代替步长内积分值。代数精度1，二阶方法。代码：
function x=tixing(a,b,n)
f=@(x) x/(4+x*x);
x=0;
h=(b-a)/n;
for k=1:n-1
x=x+f(a+k*h);
end
x=h*(f(a)+x*2+f(b))/2;
end
View Code
（复合）辛普森方法
用二次插值函数多项式近似。代数精度3，四阶方法。代码：
function t=simpson(a,b,n)
x=0;
h=(b-a)/n;
f=@(x) x/(4+x*x);
for k=0:n-1
x=x+f(a+k*h+h/2);
end
y=0;
for k=1:n-1
y=y+f(a+k*h);
end
t=h*(f(a)+4*x+2*y+f(b))/6;
View Code
（复合）科茨公式
用四次插值函数多项式近似。代数精度5，六阶方法。
龙贝格算法
近似值精度与步长选择有关，计算过程中可以通过判断当前误差调整步长以提高精度，但其计算缓慢。龙贝格算法以不同的权重组合不同区间长度的求积公式（几何意义上可以理解为在区间合适的位置划线，例如在中间画直线替代在顶上，使多出来的和少掉的抵消），从而抵消主要误差。视其误差正负号来判断组合的时候的符号。
代码：
function x=Romberg(a,b,eps)
%龙贝格算法求积分
%a,b分别为上下限，eps为精度
f=@(x) sqrt(400*sin(x)*sin(x)+100*cos(x)*cos(x));
k0=4;
k1=10;
T=zeros(4,1);
t=zeros(4,1);
T(1)=(b-a)*(f(a)+f(b))/2;
k=1;
while k<=k1
t(1)=T(1)/2;
temp=0;
for i=1:2^(k-1)
temp=temp+f(a+(2*i-1)*(b-a)/2^k);
end
temp=temp*(b-a)/2^k;
t(1)=t(1)+temp;
for i=2:min(4,k+1)
t(i)=((4^i)*t(i-1)-T(i-1))/(4^i-1);
end
if k>=k0
if abs(t(4)-T(4))<eps
x=t(4);
break;
end
end
T=t;
k=k+1;
end
View Code
第六章 常微分方程的数值解法
由于高阶常微分方程可以化为一阶常微分方程，故本章只解一节常微分方程初值问题的数值解法。答案并不需要得到精确的函数表达式，只要能获得想要点的估计值即可。
如果截断误差为步长的p+1次方，那么此方法为p阶方法。
1.欧拉方法
将函数区间n等分，得到n+1个等分点上的函数值。已知每个函数值是待求函数的导数，用直线代替该区间上的待求函数则是该线段的斜率，画出一段，再连下一段，得到一堆折线，因此又称欧拉折线法。欧拉方法是一阶方法。
代码：
function Euler(a,b,n,y0)
f1=@(x,y) 2*x/3/y/y;
f2=@(x) (1+x*x)^(1/3);
h=(b-a)/n;
x=a:h:b;
y=zeros(1,n);
y(1)=y0;
for i=1:n
y(i+1)=y(i)+h*f1(x(i),y(i));
end
x2=a:h:b;
y2=zeros(1,n);
for i=1:n+1
y2(i)=f2(x2(i));
end
plot(x,y,'-o');
hold on;
plot(x2,y2,'-ro');
for i=1:n+1
y2(i)=abs(y(i)-y2(i));
end
display(y2);
View Code
2.改进欧拉方法
欧拉方法向前作折线，也就是向前差分。改进欧拉方法向后差分，得到一个隐函数，计算量变大，稳定性增强，但仍是一阶方法。隐式欧拉公式和欧拉公式的误差主要部分符号相反，因此将它们两个加起来除以2，得到梯形公式，是二阶方法。再用显式欧拉公式算个初值，联合使用，得到改进欧拉方法，其仍为二阶方法（P151
代码：
function Euler2(a,b,n,y0)
f1=@(x,y) 2*x/3/y/y;
f2=@(x) (1+x*x)^(1/3);
h=(b-a)/n;
x=a:h:b;
y=zeros(1,n);
y(1)=y0;
for i=1:n
y(i+1)=y(i)+h/2*(f1(x(i),y(i))+f1(x(i+1),y(i)+h*(f1(x(i),y(i)))));
end
x2=a:h:b;
y2=zeros(1,n);
for i=1:n+1
y2(i)=f2(x2(i));
end
plot(x,y,'-o');
hold on;
plot(x2,y2,'-ro');
for i=1:n+1
y2(i)=abs(y(i)-y2(i));
end
display(y2);
View Code
3.龙格-库塔法
欧拉方法与改进欧拉方法的延伸，由微分中值定理引出。欧拉方法中以区间左端点斜率作为区间折线斜率，欧拉改进方法以左右端点平均值作为其斜率，龙格-库塔法想要用各种各样的方法近似这个斜率的真值，也就是找微分中值定理里那个代表不知道是啥的 ζ ，比如多取几个点，求平均值。但这还不够，最终还需要找几个参数，求得加权平均值。
代码（经典四阶方法，也被称为R-K法）：
function RungeKutta(a,b,n,y0)
h=(b-a)/n;
f1=@(x,y) 2*x/3/y/y;
f2=@(x) (1+x*x)^(1/3);
h=(b-a)/n;
x=a:h:b;
y=zeros(1,n);
y(1)=y0;
for i=1:n
xx=x(i);
yy=y(i);
k1=f1(xx,yy);
k2=f1(xx+h/2,yy+h/2*k1);
k3=f1(xx+h/2,yy+h/2*k2);
k4=f1(xx+h,yy+h*k3);
y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end
x2=a:h:b;
y2=zeros(1,n);
for i=1:n+1
y2(i)=f2(x2(i));
end
plot(x,y,'-o');
hold on;
plot(x2,y2,'-ro');
for i=1:n+1
y2(i)=abs(y(i)-y2(i));
end
display(y2);
View Code
标签：xi,end,函数,插值,存档,x2,y2,计算方法 来源： https://www.cnblogs.com/capterlliar/p/16295674.html