T1 · Lagrange 插值
依赖同目录下的 lag_.m(见本节第二段代码)。
clear; clc; close all;
rand('seed',2453655); % 学号 2453655
a = 1 + 3*rand(1,5);
b = 1 + 5*rand(1,5);
xq = 2 + 3*rand(1,2);
% --- Lagrange 插值 ---
yq = lag_(a, b, xq);
% --- 绘图 ---
xx = linspace(min(a), max(xq));
yy = lag_(a, b, xx);
plot(xx, yy, 'k-');hold on;
plot(a, b, 'ro');
plot(xq, yq, 'b*');
% Lagrange插值函数
% function yq=lag_interp(x,y,xq)
% n=length(x);
% yq=zeros(size(xq));
%
% for i=1:n
% L=ones(size(xq))
% for j=1:n
% if j~=i
% L = L .* (xq-x(j))/(x(i)-x(j))
% end
% end
% yq=yq+y(i)*L
% end
% end
function yq=lag_interp(x,y,xq)
n=length(x);
yq=zeros(size(xq));
for i=1:n
L=ones(size(xq))
for j=1:n
if j~=i
L=L.*(xq-x(j))/(x(i)-x(j))
end
end
yq=yq+y(i)*L
end
endfunction yq=lag_(x,y,xq)
n=length(x);
yq=zeros(size(xq));
for i=1:n
L=ones(size(xq));
for j=1:n
if j~=i
L=L.*(xq-x(j))/(x(i)-x(j))
end
end
yq=yq+y(i)*L;
end
endT2 · 等距节点与 Chebyshev 节点插值对比
clear; clc; close all;
rand('seed',2453655); % 学号 2453655
a = randi([5,20],1); % 区间半宽
function y=f(x)
y=1 ./(1+x .^2);
end
xx = linspace(-a, a);
c = [2 4 6 8]; % 等距节点
for k = 1:length(c)
n = c(k);
xd = linspace(-a, a, n+1); % 创建等距向量
yd = f(xd);
Lx = lagrange_interp(xd, yd, xx);
plot(xx, f(xx), 'k-'); hold on;
plot(xx, Lx, 'b-');
plot(xd, yd, 'ro');
end
% ---------- (2) Chebyshev 节点对比 ----------
n = 16; % 取较高次数凸显振荡
xd_eq = linspace(-a, a, n+1); % 等距节点
xd_ch = -a*cos((0:n)*pi/n); % Chebyshev(Lobatto) 节点
L_eq = lagrange_interp(xd_eq, f(xd_eq), xx);
L_ch = lagrange_interp(xd_ch, f(xd_ch), xx);
plot(xx, f(xx), 'k-'); hold on;
plot(xx, L_eq, 'b-');
plot(xd_eq, f(xd_eq), 'ro');
plot(xx, f(xx), 'k-'); hold on;
plot(xx, L_ch, 'b-');
plot(xd_ch, f(xd_ch), 'ro');
ylim([-0.1 1.1]);
% ===================== 局部函数 =====================
function y = lagrange_interp(xd, yd, x)
n = numel(xd);
y = zeros(size(x));
for i = 1:n
L = ones(size(x));
for j = 1:n
if j ~= i
L = L .* (x - xd(j)) / (xd(i) - xd(j));
end
end
y = y + yd(i) * L;
end
endT3 · Newton 插值
clear;clc;
rand("seed",2453655);
x=1.2:0.1:1.6;
y=1+rand(1,5);
x_dense=linspace(1.2,1.6,200);
y_dense=newton_(x,y,x_dense);
yq=newton_(x,y,1.56);
plot(x_dense,y_dense,'k-');hold on;
plot(1.56,yq,'r*');
function [yq,div_diff,table]=newton_(x,y,xq)
x=x(:);
y=y(:);
n=length(x);
table=zeros(n);
table(:,1)=y;
for j=2:n
for i=1:(n-j+1)
table(i,j)=(table(i+1,j-1)-table(i,j-1))/(x(i+j-1)-x(i))
end
end
div_diff=table(1,:);
yq=div_diff(n);
for k=(n-1):-1:1
yq=div_diff(k)+(xq-x(k)).*yq
end
endT4 · 地貌图与等高线(二维插值)
clear; clc; close all;
x = 1200:400:4000;
y = 1200:400:3600;
Z = [ ...
1130 1250 1280 1230 1040 900 500 700
1320 1450 1420 1400 1300 700 900 850
1390 1500 1500 1400 900 1100 1060 950
1500 1200 1100 1350 1450 1200 1150 1010
1500 1200 1100 1550 1600 1550 1380 1070
1500 1550 1600 1550 1600 1600 1600 1550
1480 1500 1550 1510 1430 1300 1200 980];
[Xi, Yi] = meshgrid(linspace(1200,4000,120), linspace(1200,3600,120));
[Xg, Yg] = meshgrid(x, y); % 原始观测站点坐标（与 Z 同尺寸）
methods = {'linear','nearest','spline','cubic'};
figure('Color','w','Position',[50 50 1150 780]);
for k = 1:4
Zi = interp2(x, y, Z, Xi, Yi, methods{k});
subplot(2,2,k);
surf(Xi, Yi, Zi); shading interp; hold on;
plot3(Xg, Yg, Z, 'r.', 'MarkerSize',10);
grid on;
title(['地貌图 - ', methods{k}]);
xlabel('x'); ylabel('y'); zlabel('高程 (m)');
zlim([400 1700]);
end
for k = 1:4
Zi = interp2(x, y, Z, Xi, Yi, methods{k});
subplot(2,2,k);
contourf(Xi, Yi, Zi, 15); hold on;
plot(Xg, Yg, 'r+', 'MarkerSize',6, 'LineWidth',1);
colorbar;
title(['等高线 - ', methods{k}]);
xlabel('x'); ylabel('y');
endT5 · 多项式拟合与用水量积分
clear; clc; close all;
rand('seed',2453655);
t = [0.46, 1.38, 2.40, 3.41, 4.42, 5.44,...
6.45, 7.47, 8.45, 11.49, 12.49, ...
13.42, 14.43, 15.44, 16.37, 17.38, ...
18.48, 19.50, 20.40, 23.42, 24.43, 25.45];
v = 8 + 12 * rand(22,1);
p = polyfit(t', v, 9); % 最小二乘多项式拟合
P_int = polyint(p); % 不定积分
total_water = diff(polyval(P_int, [0, 24])); % 差分求定积分
t_dense = linspace(0, 25.45, 200);
plot(t, v, 'ro', t_dense, polyval(p, t_dense), 'b-', 'LineWidth', 1.5);T6 · 多项式求根
clear; clc; close all;
rand('seed',2453655);
a = 8 + 5*rand(1);
b = 23 + rand(1);
c = randi(25,1,1);
p = [1 -a b -c];
r = roots(p);
disp(r);
% f= @(x) x .^3 - a*x .^2 + b*x - c;
% xx=linspace(min(r)-2,max(r)+2);
% plot(xx,f(xx),'k-',r,f(r),'r*');T7 · Newton 迭代法解非线性方程
clear; clc; close all;
rand('seed',2453655);
a = 1 + 5 * rand(1);
b = 0.8 + 2 * rand(1);
f = @(x) a * x.^b - exp(x);
df = @(x) a * b * x.^(b - 1) - exp(x);
% Newton 迭代法求解
x = 1.0; % 初始迭代初值 (需保证 x > 0)
tol = 1e-6; % 容差
max_iter = 100; % 最大迭代步数
for k = 1:max_iter
x_next = x - f(x) / df(x);
if abs(x_next - x) < tol
break;
end
x = x_next;
end
% xx=linspace(-3,3,200)
% plot(xx,f(xx),'b-',x,0,'r*');T8 · 数值积分(integral / trapz / 中点法)
clear;clc;close all;
rand("seed",2453655);
a=randi(10,1,1);
b=randi([5,15],1,1);
f=@(x) exp(-0.5.^x).*(x+pi/b);
n=1000;
t=linspace(0,a*pi,n+1);
w=a*pi/n;
F1=integral(f,0,a*pi);
F2=trapz(t,f(t));
F3=sum(f(t(1:n)+w/2)*w);T9 · Romberg 求积
clear;clc;close all;
a = 2.5 * rand(1, 1);
b = 3 * rand(1, 1);
f = @(x) x.^b;
tol = 1e-6; % 容差
max_k = 15; % 最大迭代层数
R = zeros(max_k, max_k);
h = a;
% 初始梯形值 R(1,1)
R(1, 1) = 0.5 * h * (f(0) + f(a));
for k = 2:max_k
nodes = (1:2^(k-2)) - 0.5;
R(k, 1) = 0.5 * R(k-1, 1) + 0.5 * h * sum(f(nodes * h));
h = h / 2;
for j = 2:k
R(k, j) = (4^(j-1) * R(k, j-1) - R(k-1, j-1)) / (4^(j-1) - 1);
end
if abs(R(k, k) - R(k-1, k-1)) < tol
break;
end
end
% 输出结果与理论值校验
I_romberg = R(k, k);
I_exact = (a^(b + 1)) / (b + 1);
fprintf('Romberg求积结果: %.8f (外推阶数: %d)\n', I_romberg, k);T10 · Jacobi 与 Gauss-Seidel 迭代
clear;clc;close all;
rand("seed",2453655);
A = [ 6, 2, 1, -2; 2, 5, 0, -2; -2, 0, 8, 5; 1, 3, 2, 7];
b = [randi(8); randi(10); -1; 0];
D = diag(diag(A)); % 先提取对角向量，再构造对角阵
L_off = -(tril(A, -1));
U_off = -(triu(A, 1));
% 构造 Jacobi 与 Gauss-Seidel 迭代矩阵并判断收敛性 (谱半径 < 1)
B_J = D \ (L_off + U_off);
B_GS = (D - L_off) \ U_off;
% 迭代求解 (统一迭代 25 步)
x_J = zeros(4,1); x_GS = zeros(4,1);
for k = 1:25
x_J = B_J * x_J + (D \ b);
x_GS = B_GS * x_GS + ((D - L_off) \ b);
end
% 输出结果对比
disp('Jacobi 解:'); disp(x_J');
disp('Gauss-Seidel 解:'); disp(x_GS');T11 · fsolve 解非线性方程组
clear;clc;close all;
F=@(t) [0.1*exp(t(1))+4*t(1)-t(2)-1;t(1)-0.125*t(1)^2-4*t(2)];
t0=[0;0];
x=fsolve(F,t0);
disp(x);T12 · 显式 Euler 与改进 Euler 法
clear;clc;close all;
a = 1 + 2 * rand(1);
t0 = 1 + rand(1);
te = 2 + rand(1);
x0 = -2;
f = @(t, x) (a * x.^2 + x) ./ t;
h_list = [0.05, 0.01];
for i = 1:length(h_list)
h = h_list(i);
t = t0:h:te;
N = length(t);
% 初始化
x_euler = zeros(1, N); x_euler(1) = x0;
x_imp = zeros(1, N); x_imp(1) = x0;
% 显式 Euler 与改进 Euler 迭代
for k = 1:N-1
% 显式欧拉法
x_euler(k+1) = x_euler(k) + h * f(t(k), x_euler(k));
% 改进欧拉法
x_pred = x_imp(k) + h * f(t(k), x_imp(k));
x_imp(k+1) = x_imp(k) + 0.5 * h * (f(t(k), x_imp(k)) + f(t(k+1), x_pred));
end
% 绘图比较
plot(t, x_euler);hold on;
plot(t, x_imp);
endT13 · ode45(Runge-Kutta)解常微分方程
clear;clc;close all;
rand("seed",2453655);
a=randi(7,1,1);
tspan=[0 0.4];
[x,y] = ode45(@(x,y) -a*y+x^2+a*x, tspan, 1);
plot(x,y,'b-',0,1,'ro');