数值分析作业

实验五 最小二乘拟合预测人口数量

x=[ 1949 1950 1952 1953 1955 1956 1957 1958 1959 1960 1961 1962 1963 1965 1966 1967 1968 1970 1971 1972 1974 1975 1976 1977 1979 1980 1982 1983 1984 ];
y=[ 5.4167 5.5196 5.7428 5.8796 6.1465 6.2828 6.4653 6.5994 6.7209 6.6207 6.5859 6.7297 6.9172 7.2538 7.4542 7.6368 7.8534 8.2992 8.5229 8.7177 9.0859 9.2420 9.3717 9.4974 9.7542 9.8705 10.1541 10.2495 10.3475 ];
xlim([1948 2000])
ylim([5 11])
plot(x,y,'*')
xlabel '年份'
ylabel '人口数（亿）'
title '散点图'
hold on
m=28;n=2; % n=3
X=zeros(n+1);
for j=1:n+1
     for i=1:n+1
          for k=1:m+1
               X(j,i)=X(j,i)+x(k)^(j+i-2);
          end
     end
end
Y=[0,0,0]; % Y=[0,0,0,0]
for j=1:n+1 
     for i=1:m+1
          Y(j)=Y(j)+y(i)*x(i)^(j-1);
     end
end

a=X\Y';
x1=[1948:1:1984]; % 画图
z=a(1)+a(2)*x1+a(3)*x1.^2; % z=a(1)+a(2)*x1+a(3)*x1.^2+a(4)*x1.^3
plot(x1,z)
legend('离散点','y=a(1)+a(2)*x+a(3)*x.^2')
title('拟合图')
 
% 预测
T=[1969 1995 2000];
for s = 1:3
    z=a(1)+a(2)*T(s)+a(3)*T(s).^2
end
% 偏差平方和Q
Q = 0;
for m = 1:29
    n=a(1)+a(2)*x(m)+a(3)*x(m).^2;
    Q = Q + (y(m)-n)^2;
end

实验六 几种数值积分方法的比较

求积函数
function X = fun(x)
    % R=6371;
    % h=439;
    % H=2384;
    % a=(2*R+H+h)/2  % 7.7825
    % c=(H-h)/2      % 972.5
    % c/a            % 0.125
    % S = 4*a*X; 
    X = sqrt(1-(0.1250*sin(x))^2);
end

复化辛普生
function S = ComplexSimpson(a, b)
% f1--四阶导函数； f2--在[a，b]范围里最大的f1
% 根据求积余项和误差限确定 h 和 n 的值
syms H x;
f1 = diff(fun(x) , 4);
f2 = max (subs(f1,x,[a:0.01:b]));
cond = abs((H/2)^4*f2*(b-a)/180) == 10E-6;
H = solve(cond, H);
h = eval(H);
n = (b-a)/h;
 

% 复化辛普生法
x1=a;
x2=a+h/2;
S = fun(a)+fun(b);
S1=0;
S2=fun(x2);
 
for k = 1:n-1
    x1=x1+h;
    x2=x2+h;
    S1=S1+fun(x1);
    S2=S2+fun(x2);
end 
S = h*(S+2*S1+4*S2)/6;
S = 4*7.7825*S;
end


龙贝格
function S = Romberg(a,b,e)
% 龙贝格
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a; 
T(1,1)=h*(fun(a)+fun(b))/2;
d=b-a; %误差间距
while e<=d
    k=k+1;
    h=h/2;
    sum=0;
    % 计算第一列T
    for i=1:n
        sum=sum+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*sum;
    % 迭代
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);
    end
    n=n*2;
    d=abs(T(k+1,k+1)-T(k,k));
end
S = T(k+1,k+1);
S = 4*7.7825*S; 





复化两点高斯—勒让德
function S = Gauss(a,b,n)
% 输入积分上下限b、a，区间分段数n
h=(b-a)/n;
S=0;
c=1/(2*sqrt(3));
for k = 0:n-1
    S=S+fun(a+(k+1/2-c)*h)+fun(a+(k+1/2+c)*h);
end
S = h*S/2;
S = 4*7.7825*S;
end
    
调用函数
ComplexSimpson(0, pi/2)
Romberg(0,pi/2,10E-6)
Gauss(0,pi/2,3)

实验七 数值求解常微分方程初值问题

欧拉方法
x = 1;
y = 0;
y1 = 2*y/x+x^2*exp(x); % (1,2]  1+0.1*k
t = eye(40,1);
for k = 1:40
    x = 1+0.025*k;
    y = y+0.025*y1;
    y1 = 2*y/x+x^2*exp(x);
    Y = x^2*(exp(x)-exp(1));
    t(k) = Y-y;
end

for k = 1:10
    x = 1+0.1*k;
    y = y+0.1*y1
    y1 = 2*y/x+x^2*exp(x);
    Y = x^2*(exp(x)-exp(1));
    t(k) = Y-y;
end


改进的欧拉方法
y = 0;
x = [1:0.1:2]; % 0.1
t = eye(20,1);
for k = 1:10
    y2 = y+0.1*(2*y/x(k)+x(k)^2*exp(x(k))); % 0.1
    y3 = y+0.1*(2*y2/x(k+1)+x(k+1)^2*exp(x(k+1))); % 0.1
    y = (y2+y3)/2;
    Y = x(k+1)^2*(exp(x(k+1))-exp(1));
    t(k) = Y-y;
End

四阶标准龙格—库塔法
h = 0.1;
x = [1:0.1:2];
y = eye(11,1);
y(1) = 0;
t = eye(10,1);
 
for k = 1:10
    k1 = 2*y(k)/x(k)+x(k)^2*exp(x(k));
    k2 = 2*(y(k)+h*k1/2)/(x(k)+h/2)+(x(k)+h/2)^2*exp(x(k)+h/2);
    k3 = 2*(y(k)+h*k2/2)/(x(k)+h/2)+(x(k)+h/2)^2*exp(x(k)+h/2);
    k4 = 2*(y(k)+h*k3)/(x(k)+h)+(x(k)+h)^2*exp(x(k)+h);
    y(k+1) = y(k) + h*(k1+2*k2+2*k3+k4)/6;
    Y = x(k+1)^2*(exp(x(k+1))-exp(1));
    t(k+1) = Y-y(k+1);
end