数值分析作业

问题：

分别使用Lagrange和Newton插值公式计算函数的值。

对于Newton插值要求先输出差商表，在计算x=7.5时4次插值多项式的值，并与四次Lagrange插值结果进行比较

x	4	5	6	7	8
y	2.0000	2.23607	2.44949	2.64575	2.82843

Lagrange

function y = Lagrange(a,b,x,n)
%Lagrange插值法 a,b分别为x以及对应的f(x),x为插入数据，n为插值多项式次数
% a=4:8
% b=[2 2.23607 2.44949 2.64575 2.82843]
if length(a)~=length(b)
    error('输入数据不同维')
end
if n>=length(a)
    error('n值太大')
end
d=abs(x-a);
[~,idx]=sort(d);
S=0;
for i=1:n+1
    t1=1;
    t2=1;
    for j=1:n+1
        if j~=i
            t1=t1*(x-a(idx(j)));
            t2=t2*(a(idx(i))-a(idx(j)));
        end
    end
    S=S+t1/t2*b(idx(i));
end
y=S;
end

Newton

•  Chashang_List

function A = Chashang_List(a,b)
%Newton插商计算函数
if length(a)~=length(b)
    error('输入数据不同维')
end
N=length(a);
A=zeros(N,N+1);
A(:,1)=a;
A(:,2)=b;
for j=3:N+1
    for i=j-1:N
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(A(i,1)-A(i-(j-2),1));
    end
end
ChaShang_List=A
end

---

function y = Newton(a,b,x,n)
% Newton插值法 a,b分别为x以及对应的f(x),x为插入数据，n为插值多项式次数
% a=4:8
% b=[2 2.23607 2.44949 2.64575 2.82843]
if n>=length(a)
    error('n值太大')
end
A=Chashang_List(a,b);
S=A(1,2);
X=cumprod(x-a);
for i=1:n
    S=S+A(i+1,i+2)*X(i);
end
y=S;
end

欧拉 

function y = Euler (fun, x0, y0, h, b)
n = (b - x0)/h;
A = zeros (n+1, 4);
A(1,1) = x0;
A(1,2) = y0;
A(1,1) = y0;
for i = 1 : n 
A(i+1,1) = A (i,1) + h;
A(i+1,2) = A (i,2) + h*fun (A(i,1), A(i,2));
A(i+1,3) = A (i,2) + h*fun (A(i+1,1), A(i+1,2));
yp = A(i,4) + h*fun (A (i,1), A(i,4));
yc = A(i,4) + h*fun (A (i+1,1), yp);
A (i+1,4) = 1/2*(yp + yc);
end 
y = A;


fun=@(x,y)-y+x+1
y = Euler(fun, 0, 1, 0.5, 1)

y =

    1.0000    1.0000         0         0
    1.5000    1.5000    1.5000    0.8750
    2.0000    2.0000    2.0000    1.6094

龙贝格

function y = Runge_Kutta (fun, x0, y0, h, b)
n = (b-x0)/h;
A = zeros (n+1, 6);
A(1,1) = x0;
A(1,6) = y0;

for i = 1 : n 
A(i+1,1) = A (i,1) + h;  
A(i+1,2) = fun (A(i,1), A(i,2));  %k1
A(i+1,3)=fun (A(i,1)+h/2, A(i, 6)+h/2* A(i+1,2));  %k2
A(i+1,4)=fun (A(i,1)+h/2, A(i, 6)+h/2* A(i+1,3));  %k3
A(i+1,5)=fm(A(i,1)+h, A(i,6)+h* A(i+1,4));  %k4
A(i+1,6)=A(i,6)+h/6*(A(i+1,2)+2*A(i+1,3)+2*A(i+1,4)+A(i+1,5));
end 
y = A;


fun=@(x,y)-3*y+2*x+2
y = Euler(fun, 0, 1, 0.1, 1)

y =

    1.0000    1.0000         0         0
    1.1000    1.1000    1.0900    0.3500
    1.2000    1.1900    1.1830    0.6278
    1.3000    1.2730    1.2681    0.8517
    1.4000    1.3511    1.3477    1.0355
    1.5000    1.4258    1.4234    1.1894
    1.6000    1.4980    1.4964    1.3211
    1.7000    1.5686    1.5675    1.4362
    1.8000    1.6380    1.6372    1.5390
    1.9000    1.7066    1.7061    1.6326
    2.0000    1.7746    1.7742    1.7193

 用迭代法求方程 e^x +10x - 2 = 0 的近似根， 取初值 x = 0 ， 迭代公式，要求精度为0.5 *10^-3 ;

X0=0;
errf=1;
i=0;
while errf>0.5*10^(-3)
i=i+1;
x=(2-exp(X0))/10;
errf=abs((x-X0)/X0);
x0=x;
end

插值

function yy = mymulNewtonCotes(ft,a,b,m,n)

%复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和,

%参考的输入形式为 mymulNewtonCotes(ft,0,1,10,2)

%参数说明:

%    ft 被积函数，此题中ft=@(t)t.*exp(t^2/2)

%    a 积分下限

%    b 积分上限

%    m 将区间[a,b]等分的子区间数量

%    n 采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性

%     (1)n=1时为复化梯形公式

%     (2)n=2时为复化辛普森公式

function yy = mymulNewtonCotes(ft,a,b,m,n)

XX= linspace(a,b,m+1);
for l= 1:m
    s(1)= myNewtonCotes(ft,xx(1),xx(1+1),n);
end
yy = sum(s);

function [y,Ck,Ak]= myNewtonCotes(ft,a,b,n)%牛顿-科特斯数值积分公式

% Ck——科特斯系数
% Ak———求积系数
% y—牛顿-科特斯数值积分结果

xk = linspace(a,b,n+1);
for j = 1:n+1
ff(j)= ft(xk(j));
end

%计算科特斯系数
for i=1:n+1
    k=i-1;
    Ck(i)=(-1)y(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n);
end

%计算求积系数
Ak=(b-a)*Ck;

%求和算积分
y=Ak*ff;

function f=intfun(t,n,k)
%科特斯系数中的积分表达式

f=1;
for i=[0:k-1,k+1:n]
    f=f.*(t-i);
end