matlab编程练习

x = [0.2, 0.4, 0.6, 0.8, 1.0];
y = [0.98, 0.92, 0.81, 0.64, 0.38];
n = length(x);
F = zeros(n, n);
F(:,1) = y';

for j = 2:n
    for i = j:n
        F(i,j) = (F(i,j-1) - F(i-1,j-1)) / (x(i) - x(i-j+1));
    end
end

newton_interp = @(xx) F(1,1) + ...
    (xx - x(1))*(F(2,2) + ...
    (xx - x(2))*(F(3,3) + ...
    (xx - x(3))*(F(4,4) + ...
    (xx - x(4))*F(5,5))));
pp = csape(x, y, 'variational');
x_new = 0.2:0.08:1.0;
y_newton = arrayfun(newton_interp, x_new);
y_spline = fnval(pp, x_new);

plot(x, y, 'o', x_new, y_newton, '-', x_new, y_spline, '--');
legend('原始数据', '牛顿插值', '三次样条');
title('插值结果对比');
xlabel('x'); ylabel('f(x)');

t = 0:24;
x = [15,14,14,14,14,15,16,18,20,22,23,25,28,31,32,31,29,27,25,24,22,20,18,17,16];
po = polyfit(t, x, 2);
y2 = polyval(po, t);
error2 = sum((x - y2).^2);
py = polyfit(t, x, 3);
y3 = polyval(py, t);
error3 = sum((x - y3).^2);
pm = polyfit(t, x, 4);
y4 = polyval(pm, t);
error4 = sum((x - y4).^2);
model = @(params, t) params(1)*exp(-params(2)*(t - params(3)).^2);
params0 = [32, 0.1, 14];
[params, resnorm] = lsqcurvefit(model, params0, t, x);
y_gauss = model(params, t);
error_gauss = resnorm;
figure;
plot(t, x, 'ko', 'DisplayName', '原始数据');
hold on;
plot(t, y2, 'b-', 'DisplayName', sprintf('二次 (误差=%.2f)', error2));
plot(t, y3, 'r-', 'DisplayName', sprintf('三次 (误差=%.2f)', error3));
plot(t, y4, 'g-', 'DisplayName', sprintf('四次 (误差=%.2f)', error4));
plot(t, y_gauss, 'm--', 'DisplayName', sprintf('高斯 (误差=%.2f)', error_gauss));
legend;
title('气温变化拟合比较');
xlabel('时间 t'); ylabel('温度 x(t)');

a = 0; b = pi/2;
simp = (b - a)/6 * (sin(a) + 4*sin((a + b)/2) + sin(b));
fprintf('辛普森公式结果: %.8f\n', simp);
n = 4;
h = (b - a)/n;
x = a:h:b;
sum_odd = sum(sin(x(2:2:end-1)));
sum_even = sum(sin(x(3:2:end-2)));
comp_simp = h/3 * (sin(a) + sin(b) + 4*sum_odd + 2*sum_even);
fprintf('复合辛普森结果(n=4): %.8f\n', comp_simp);

A = [3, -1, 2, 7;
     -1, 2, -2, -1;
      2, -2, 4, 0];
n = size(A,1);
for k = 1:n-1
    for i = k+1:n
        m = A(i,k)/A(k,k);
        A(i,k:n+1) = A(i,k:n+1) - m*A(k,k:n+1);
    end
end
x = zeros(n,1);
x(n) = A(n,n+1)/A(n,n);
for i = n-1:-1:1
    x(i) = (A(i,n+1) - A(i,i+1:n)*x(i+1:n)) / A(i,i);
end

fprintf('解: x1=%.6f, x2=%.6f, x3=%.6f\n', x(1), x(2), x(3));

A = [8, -3, 2;
     4, 11, -1;
     6, 3, 12];
b = [20; 33; 36];
x = A\b;
fprintf('解: x1=%f, x2=%f, x3=%f\n', x(1), x(2), x(3));