超级会员免费看
MATLAB 数值计算与线性代数问题求解
1. 多项式拟合与插值问题
1.1 直线拟合
输入四个值代表点 $(x_1, y_1)$ 和 $(x_2, y_2)$，将其分别存储在向量 $x$ 和 $y$ 中，$x$ 存储 $x$ 坐标，$y$ 存储 $y$ 坐标。使用
polyfit
命令拟合直线 $y = mx + c$，其中梯度 $m$ 是 $p$ 的第一个元素，截距 $c$ 是第二个元素。最后显示结果。若用户输入的两个 $y$ 值相同，则 $m = 0$，$c$ 等于该值；若两个 $x$ 值相同，MATLAB 会尝试给出无穷梯度的直线，形式为 $x = d$。
% 示例代码
% 假设输入的四个值
x1 = 1; y1 = 2; x2 = 3; y2 = 4;
x = [x1, x2];
y = [y1, y2];
p = polyfit(x, y, 1);
m = p(1);
c = p(2);
disp(['直线方程为 y = ', num2str(m), 'x + ', num2str(c)]);
1.2 二次多项式拟合
二次多项式形式为 $y = a + b(x - x_0) + c(x - x_0)(x - x_1)$，选择一个点的 $x$ 坐标作为 $x_0$，为方便起见，取原点为第一个点，即 $x_0 = 0$，此时 $a = 0$。通过其他点确定 $b$ 和 $c$ 的值。例如，使用点 $(2, -1)$ 可得 $-1 = 2b$，再结合另一个点可求出 $c$ 的值。
% 示例代码
% 已知条件
x0 = 0;
% 使用点 (2, -1)
x1 = 2; y1 = -1;
b = -1/2;
% 假设另一个点确定 c 的值
% 这里省略具体计算过程
c = 1/2;
syms x;
y = -1/2*x + 1/2*x*(x - 2);
disp(['二次多项式方程为 y = ', char(y)]);
1.3 插值误差计算
使用以下代码计算插值误差：
x = 0:10;
co = [1 3 2];
y = x.^2+3*x+2;
for i = 1:3
xv = i-0.5;
p = polyfit(x(i:(i+1)),y(i:(i+1)),1);
err(i) = polyval(p,xv)-polyval(co,xv);
end
disp(['插值误差为 ', num2str(err)]);
该代码在每种情况下给出的误差均为 $1/4$，这可通过理解插值方法的误差（考虑二阶导数）来预期。
1.4 特定二次多项式
二次多项式形式为 $cx(\pi - x)$（因为在 $0$ 和 $\pi$ 处为零），通过曲线经过的最后一个点确定 $c$ 的值，得到 $y = \frac{4x}{\pi^2}(\pi - x)$。
1.5 三次多项式
三次多项式在 $0$ 和 $\pi$ 处为零，形式为 $x(\pi - x)(a(x - \pi/2) + b)$，通过其他点确定 $a$ 和 $b$ 的值。例如，点 $(\pi/2, 1)$ 给出 $b = \frac{4}{\pi^2}$，点 $(-\pi/2, -1)$ 给出 $a = \frac{16}{3\pi^2}$，最终得到三次多项式。
% 示例代码
syms x pi;
a = 16/(3*pi^3);
b = 4/pi^2;
y = x*(pi - x)*(a*(x - pi/2) + b);
disp(['三次多项式方程为 y = ', char(y)]);
1.6 样条插值
使用 MATLAB 的
spline
命令进行样条插值：
x = -pi:(pi/2):pi;
y = [0 -1 0 1 0];
z = -pi:(pi/20):pi;
f = spline(x,y,z);
plot(z,f,x,y,'o','MarkerSize',14)
1.7 最小二乘法拟合
使用最小二乘法拟合函数 $a\sin x + b\cos x$ 到给定数据。首先计算误差平方和：
$e = \sum_{i = 1}^{n}(a\sin x_i + b\cos x_i - f_i)^2$
对 $a$ 和 $b$ 求偏导数：
$\frac{\partial e}{\partial a} = \sum_{i = 1}^{n}\sin x_i(a\sin x_i + b\cos x_i - f_i)$
$\frac{\partial e}{\partial b} = \sum_{i = 1}^{n}\cos x_i(a\sin x_i + b\cos x_i - f_i)$
将这些方程重写为矩阵形式：
$\begin{pmatrix}
\sum_{i = 1}^{n}\sin^2 x_i & \sum_{i = 1}^{n}\sin x_i\cos x_i \
\sum_{i = 1}^{n}\cos x_i\sin x_i & \sum_{i = 1}^{n}\cos^2 x_i
\end{pmatrix}
\begin{pmatrix}
a \
b
\end{pmatrix}
=
\begin{pmatrix}
\sum_{i = 1}^{n}f_i\sin x_i \
\sum_{i = 1}^{n}f_i\cos x_i
\end{pmatrix}$
使用以下代码求解：
x = 0:0.1:1.0;
f = [3.16 3.01 2.73 2.47 2.13 1.82 1.52 1.21 0.76 0.43 0.03];
A = [sum(sin(x).^2) sum(cos(x).*sin(x)); sum(cos(x).*sin(x)) sum(cos(x).^2)];
r = [sum(f.*sin(x)); sum(f.*cos(x))];
sol = A\r;
disp(['最小二乘法拟合结果 a = ', num2str(sol(1)), ', b = ', num2str(sol(2))]);
1.8 样条插值系数与误差
使用以下代码进行样条插值并计算误差：
x = -pi:(pi/2):pi;
y = [0 -1 0 1 0];
z = -pi:(pi/10):pi;
pp = spline(x,y);
f = spline(x,y,z);
plot(z,f,x,y,'o','MarkerSize',14)
true = sin(z);
err = sum((true-f).^2);
disp(['样条插值误差为 ', num2str(err)]);
disp(['样条插值系数为 ']);
disp(pp.coefs);
1.9 修正代码
x = 2:11;
f = polyval([1 0 0 -1],x) + sin(x);
% x = 4.5
r = 3:4;
c = polyfit(x(r),f(r),1);
yy = polyval(c,4.5);
disp(['x = 4.5 时的预测值为 ', num2str(yy)]);
% x = 15 (extrapolation)
r = length(x)-1:length(x);
c = polyfit(x(r),f(r),1);
yy = polyval(c,15);
disp(['x = 15 时的预测值为 ', num2str(yy)]);
2. 线性代数问题
2.1 矩阵运算
使用以下代码进行矩阵乘法、加法、转置等运算：
A = [3 0 -1; -4 2 2];
B = [-1 7; 3 5; -2 0];
C = [2 0; -1 -3];
disp('A * B = ');
disp(A*B);
disp('B * A = ');
disp(B*A);
disp('A + transpose(B) = ');
disp(A+transpose(B));
disp('A * C = ');
disp(A*C);
disp('A * transpose(C) = ');
disp(A*transpose(C));
disp('3 * C + 2 * transpose(A*B) = ');
disp(3*C+2*transpose(A*B));
disp('(A * B) * C = ');
disp((A*B)*C);
disp('A * (B * C) = ');
disp(A*(B*C));
2.2 矩阵元素设置
r = 1:4;
A = zeros(4);
A(1,r) = r;
A(r,4) = flipud(r');
disp('矩阵 A = ');
disp(A);
2.3 特殊矩阵创建
a = diag(ones(1,9),1)+diag(-ones(1,9),-1);
disp('特殊矩阵 a = ');
disp(a);
2.4 矩阵运算示例
A = [ 1 2; 3 4];
B = [3 4; -1 2];
disp('A * B = ');
disp(A*B);
C = [3 5; 6 -2];
D = [-1 0; 2 1];
disp('2 * C - 4 * D = ');
disp(2*C-4*D);
E = [1 3 5];
F = [2 -1; -1 0; 7 -2];
disp('E * F = ');
disp(E*F);
2.5 矩阵交换性
考虑一般矩阵 $A = \begin{pmatrix}
a & b \
c & d
\end{pmatrix}$ 与矩阵 $X = \begin{pmatrix}
\alpha & \beta \
\beta & \alpha
\end{pmatrix}$ 相乘，通过比较 $XA$ 和 $AX$ 可知，当 $\beta \neq 0$ 时，$a = d$ 且 $b = c$ 时两矩阵可交换。以下代码可验证：
stl = 'Top left element of matrix ';
sbl = 'Bottom left element of matrix ';
for j = 1:2
a(j) = input([stl num2str(j), ': ']);
b(j) = input([sbl num2str(j), ': ']);
end
A = [a(1) b(1); b(1) a(1)];
B = [a(2) b(2); b(2) a(2)];
disp('A * B = ');
disp(A*B);
disp('B * A = ');
disp(B*A);
2.6 对称与反对称矩阵
设矩阵 $B$ 的 $(i, j)$ 元素为 $a_{i,j} + a_{j,i}$，则 $B$ 是对称矩阵；矩阵 $C$ 的 $(i, j)$ 元素为 $a_{i,j} - a_{j,i}$，则 $C$ 是反对称矩阵。
2.7 旋转矩阵
旋转矩阵可通过以下代码构建：
theta = 0;
A0 = [cos(theta) sin(theta); -sin(theta) cos(theta)];
theta = pi/2;
A1 = [cos(theta) sin(theta); -sin(theta) cos(theta)];
theta = pi;
A2 = [cos(theta) sin(theta); -sin(theta) cos(theta)];
disp('theta = 0 时的旋转矩阵 A0 = ');
disp(A0);
disp('theta = pi/2 时的旋转矩阵 A1 = ');
disp(A1);
disp('theta = pi 时的旋转矩阵 A2 = ');
disp(A2);
旋转矩阵的作用是将点旋转 $\theta$ 弧度，其逆矩阵为 $A^{-1} = \begin{pmatrix}
\cos(-\theta) & \sin(-\theta) \
-\sin(-\theta) & \cos(-\theta)
\end{pmatrix} = \begin{pmatrix}
\cos\theta & -\sin\theta \
\sin\theta & \cos\theta
\end{pmatrix}$，可通过矩阵乘法验证。
2.8 线性方程组求解
2.8.1 二元线性方程组
A = [3 4; -1 2];
b = [2 ; 0];
x = A\b;
disp(['二元线性方程组的解为 x = ', num2str(x(1)), ', y = ', num2str(x(2))]);
2.8.2 三元线性方程组
A = [1 1 2; 1 -1 -3; -2 -5 1];
b = [1; 0 ; 4];
x = A\b;
disp(['三元线性方程组的解为 x = ', num2str(x(1)), ', y = ', num2str(x(2)), ', z = ', num2str(x(3))]);
2.9 矩阵加法与乘法
矩阵 $A$ 和 $B$ 可相加，因为它们大小相同（三行两列）；矩阵 $A$ 和 $C$ 可相乘，因为 $A$ 的列数等于 $C$ 的行数。
a = [1 -1; 0 2; 3 2];
b = [2 -1; -1 0; 3 2];
c = [-1 0; 2 1];
disp('a + b = ');
disp(a+b);
disp('a * c = ');
disp(a*c);
disp('(a - b) * c = ');
disp((a-b)*c);
disp('a * c - b * c = ');
disp(a*c-b*c);
2.10 矩阵乘法示例
disp('[1 -1 2; 3 0 1] * [3; 2; 1] = ');
disp([1 -1 2; 3 0 1]*[3; 2; 1]);
disp('[5 -2; -1 2] * [4 0 1 -1; 2 1 -2 -1] = ');
disp([5 -2; -1 2]*[4 0 1 -1; 2 1 -2 -1]);
2.11 矩阵与单位矩阵乘法
矩阵与单位矩阵相乘结果不变。
2.12 矩阵转置与乘法
A = [3 2 -1; 0 -1 -2];
disp('A * transpose(A) = ');
disp(A*transpose(A));
disp('transpose(A) * A = ');
disp(transpose(A)*A);
2.13 行向量与转置
一般行向量 $(x_1, x_2, \cdots, x_N)$ 的转置为列向量 $\begin{pmatrix}
x_1 \
x_2 \
\vdots \
x_N
\end{pmatrix}$，则 $xx^T = x_1^2 + x_2^2 + \cdots + x_N^2$，结果为正的标量。
2.14 矩阵方程
矩阵方程可展开为线性方程组，例如：
$\begin{pmatrix}
1 & 1 & 1 \
1 & -2 & -1 \
-1 & 3 & -1
\end{pmatrix}
\begin{pmatrix}
x \
y \
z
\end{pmatrix}
=
\begin{pmatrix}
0 \
2 \
-1
\end{pmatrix}$
2.15 系统特征判断
a = [3 2; 3 -2]; b=[7; 7];
% 假设 solns 函数已定义
% solns(a,b)
a = ones(6);
for r = 2:6
a(r,r) = -1;
end
b = ones(6,1);
% solns(a,b)
2.16 矩阵方程组求解
A = [1 0 0 -1; -1 2 -1 0; 0 -1 2 -1; 0 0 0 1];
r = [0 1;0 0; 0 0; 1 0];
sols = A\r;
disp('矩阵方程组的解为 ');
disp(sols);
2.17 矩阵行列式
s = pi:pi/3:(2*pi);
ns = length(s);
for j = 1:ns
ss = s(j);
A = [0 1 ss; ss 0 1; 1 ss 0];
z(j) = det(A);
end
c = polyfit(s,z,3);
disp(['矩阵行列式的多项式系数为 ', num2str(c)]);
2.18 矩阵幂运算
n = input('What power :');
b = [0 1; -1 0];
switch mod(n,4)
case 0
bn = eye(2);
case 1
bn = b;
case 2
bn = -eye(2);
case 3
bn = -b;
end
disp(['矩阵 b 的 ', num2str(n), ' 次幂为 ']);
disp(bn);
2.19 特征值计算
a = [1 0 0 -1; 0 1 0 0; 0 0 1 0; -1 0 0 1];
disp('矩阵 a 的特征值为 ');
disp(eig(a));
2.20 矩阵幂的归纳证明
通过归纳法证明 $A^n = PD^nP^{-1}$。
2.21 特征多项式与特征值验证
co = charpoly(a);
disp('矩阵 a 的特征多项式的根（即特征值）为 ');
disp(roots(co));
2.22 特定方程的特征值与解
方程的特征值为 $\frac{3 \pm \sqrt{5}}{2}$，解为 $x(t) = \frac{1}{\sqrt{5}}[\frac{1}{2}((3 + \sqrt{5})e^{\frac{(3 - \sqrt{5})t}{2}} - (3 - \sqrt{5})e^{\frac{(3 + \sqrt{5})t}{2}})I + (e^{\frac{(3 - \sqrt{5})t}{2}} - e^{\frac{(3 + \sqrt{5})t}{2}})A]\begin{pmatrix}
1 \
-1
\end{pmatrix}$。
3. 积分问题
3.1 函数值生成
for i = 1:12
switch mod(i,3)
case 0
f(i) = 1;
case 1
f(i) = 2;
case 2
f(i) = 3;
end
end
disp('函数值 f = ');
disp(f);
3.2 积分权重计算
3.2.1 三分之一规则
当 $N = 9$ 时：
-
rodd=1:2:N
得到
[1 3 5 7 9]
-
reven=2:2:(N-1)
得到
[2 4 6 8]
-
weights(rodd=2)
得到
[2 0 2 0 2 0 2 0 2]
-
weights(1)=1
得到
[1 0 2 0 2 0 2 0 2]
-
weights(N)=1
得到
[1 0 2 0 2 0 2 0 1]
-
weights(reven)=4
得到
[1 4 2 4 2 4 2 4 1]
3.2.2 八分之三规则
当 $N = 10$ 时：
-
m=(N-1)/3
得到
3
-
rdiff=3*(1:(m-1))+1
得到
[4 7]
-
weights=3*ones(1,N)
得到
[3 3 3 3 3 3 3 3 3 3]
-
weights(1)
得到
[1 3 3 3 3 3 3 3 3 3]
-
weights(N)
得到
[1 3 3 3 3 3 3 3 3 1]
-
weights(rdiff)=2
得到
[1 3 3 2 3 3 2 3 3 1]
3.3 自定义函数
function [val] = fn(x)
val = log(x+sqrt(x.^2+1));
end
3.4 积分计算
N = 40;
x = linspace(1,3,N);
f = x.^2-3*x+2;
h = x(2)-x(1);
integral = (sum(f)-f(1)/2-f(N)/2)*h;
disp(['积分结果为 ', num2str(integral)]);
3.5 辛普森三分之一规则积分
x = linspace(0,1,11);
h = x(2)-x(1);
N = length(x);
rodd = 1:2:N;
reven = 2:2:(N-1);
weights(rodd) = 2; weights(1) = 1;
weights(N) = 1; weights(reven) = 4;
f = x.^3-x+1;
integral = h/3*sum(weights.*f);
disp(['辛普森三分之一规则积分结果为 ', num2str(integral)]);
3.6 不同点数的积分误差
Ns = 5:2:19;
for N = Ns
clear rodd reven weights f x
x = linspace(0,pi,N);
h = x(2)-x(1);
rodd = 1:2:N;
reven = 2:2:(N-1);
weights(rodd) = 2;
weights(1) = 1; weights(N) = 1;
weights(reven) = 4;
f = sin(x);
integral(N) = h/3*sum(weights.*f);
end
plot(Ns,abs(integral-2))
xlabel('\theta degrees')
ylabel('误差')
3.7 截断积分
X = input('Truncate at:');
N = ceil(X)*3;
x = linspace(0,X,N);
h = x(2)-x(1);
f = 1./sqrt(x.^2+1);
int = (sum(f)-f(1)/2-f(N)/2)*h;
disp(['截断积分结果为 ', num2str(int)]);
3.8 弧长计算
theta = 0:pi/20:(pi/2-pi/20);
N = 20;
for it = 1:length(theta);
theta1 = theta(it);
clear grid f
grid = linspace(theta1,pi-theta1,N);
f = sqrt(1+cos(grid).^2);
h = grid(2)-grid(1);
arclen(it) = (sum(f)-f(1)/2-f(N)/2)*h;
end
plot(theta/pi*180,arclen)
xlabel('\theta degrees')
ylabel('Arc length')
3.9 积分分析与计算
clear all
epsil = input('Epsilon :');
int1 = 2*epsil^(0.5)-epsil^(2.5)/5;
N = 100;
x = linspace(epsil,10,N);
h = x(2)-x(1);
f = cos(x)./sqrt(x);
int2 = (sum(f)-f(1)/2-f(N)/2)*h;
int = int1+int2;
disp(['积分结果为 ', num2str(int)]);
3.10 二次多项式积分
二次多项式 $y(x) = a_0 + (x - x_0)\Delta a_1 + (x - x_0)(x - x_1)a_2$，通过积分 $\int_{x_0}^{x_2} y(x) dx$ 计算。
3.11 积分函数与计算
function [f] = fxlnx(x)
f = x.*log(x);
end
result = quad('fxlnx',1,2);
disp(['积分结果为 ', num2str(result)]);
4. 微分方程问题
4.1 一阶微分方程
4.1.1 精确解
将方程 $\frac{1}{y}\frac{dy}{dt} = -\sqrt{t}$ 两边积分得到 $\ln y = -\frac{2}{3}t^{\frac{3}{2}} + C$，再根据初始条件 $y(0) = 1$ 确定 $A = 1$，精确解为 $y = e^{-\frac{2}{3}t^{\frac{3}{2}}}$。
4.1.2 数值解
dt = 0.05;
t = 0.0:dt:1.0;
y = zeros(size(t));
y(1) = 1;
for ii=1:(length(t)-1)
y(ii+1) = y(ii) + dt * (-y(ii)*sqrt(t(ii)));
end
exact = exp(-2/3*(t).^(3/2));
plot(t,y,t,exact,'--')
4.2 一阶微分方程的数值解与误差
考虑方程 $\frac{1}{y}\frac{dy}{dt} = \pm t$，精确解为 $y = e^{\pm\frac{t^2}{2}}$。数值解的递推公式为 $y_{n + 1} = y_n(1 \pm n\Delta t^2)$。通过计算不同 $n$ 时的数值解，并与精确解比较，得到绝对误差和相对误差。
4.3 一阶微分方程数值解
dt = pi/10;
t = 0.0:dt:10.0*pi;
y = zeros(size(t));
y(1) = 0;
for ii=1:(length(t)-1)
y(ii+1) = y(ii) + dt * (sin(t(ii))+sin(y(ii)));
end
plot(t,y)
4.4 一阶线性微分方程
4.4.1 精确解
通过乘以积分因子 $e^{\frac{t}{3}}$ 并积分，结合边界条件得到精确解 $y = \frac{9}{2}(1 - e^{-\frac{t}{3}}) - \frac{3t}{2}$。
4.4.2 数值解
clear all
dt = 1/3;
t = 0.0:dt:1;
y = zeros(size(t));
y(1) = 0;
for ii = 1:(length(t)-1)
y(ii+1) = 1/(1+dt/3)*(y(ii)-dt*t(ii+1)/2);
end
ts = t; ys = y;
dt = 1/1000;
t = 0.0:dt:1;
y = zeros(size(t));
y(1) = 0;
for ii = 1:(length(t)-1)
y(ii+1) = 1/(1+dt/3)*(y(ii)-dt*t(ii+1)/2);
end
exact = 9/2*(1-exp(-t/3))-3/2*t;
plot(ts,ys,'*',t,exact,t,y)
4.5 一阶微分方程求解器
function [value] = odes(t,y)
value = t^2-y^2;
end
y0 = 0;
tspan = [0 2];
[t,y] = ode45('odes',tspan,y0);
plot(t,y)
4.6 一阶线性微分方程
4.6.1 精确解
通过乘以积分因子 $e^t$ 并积分，结合初始条件得到精确解 $y(t) = t^2 - 2t + 2 - e^{-t}$。
4.6.2 数值解
N = 20;
t = linspace(0,2,N);
dt = t(2)-t(1);
y(1) = 1;
for j = 1:(N-1)
y(j+1) = y(j)+dt*(-y(j)+t(j)^2);
end
ex = t.^2-2*t+2-exp(-t);
plot(t,y,t,ex,'*')
4.7 二阶微分方程
4.7.1 问题一
精确解为 $y(x) = 2\sin x - x\cos x + x(-2\sin 1 + \cos 1 - 1) + 1$。
x = 0.0:0.1:1.0;
N = length(x);
h = x(2)-x(1);
a = 1/h^2*ones(size(x));
b = -2/h^2*ones(size(x));
c = 1/h^2*ones(size(x));
r = x.*cos(x);
a(1) = 0; b(1) = 1; r(1) = 1;
c(N) = 0; b(N) = 1; r(N) = 0;
% Forward sweep
for j = 2:N
b(j) = b(j)-c(j)*a(j-1)/b(j-1);
r(j) = r(j)-c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
y(j) = r(j)/b(j)-a(j)*y(j+1)/b(j);
end
plot(x,y)
4.7.2 问题二
精确解为 $y(x) = 2\sin x - x\cos x - x - 2\sin 1 + \cos 1 + 1$。
x = 0.0:0.1:1.0;
N = length(x);
h = x(2)-x(1);
a = 1/h^2*ones(size(x));
b = -2/h^2*ones(size(x));
c = 1/h^2*ones(size(x));
r = x.*cos(x);
a(1) = -1; b(1) = 1; r(1) = 0;
c(N) = 0; b(N) = 1; r(N) = 0;
% Forward sweep
for j = 2:N
b(j) = b(j)-c(j)*a(j-1)/b(j-1);
r(j) = r(j)-c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
y(j) = r(j)/b(j)-a(j)*y(j+1)/b(j);
end
plot(x,y)
4.7.3 问题三
精确解为 $y(x) = \frac{e^{-2x} - 1}{e^{-2} - 1}$。
x = 0.0:0.1:1.0;
N = length(x);
h = x(2)-x(1);
a = (1/h^2+2/(2*h))*ones(size(x));
b = -2/h^2*ones(size(x));
c = (1/h^2-2/(2*h))*ones(size(x));
r = 0;
a(1) = 0; b(1) = 1; r(1) = 0;
c(N) = 0; b(N) = 1; r(N) = 1;
% Forward sweep
for j = 2:N
b(j) = b(j)-c(j)*a(j-1)/b(j-1);
r(j) = r(j)-c(j)*r(j-1)/b(j-1);
end
% Final equation
y(N) = r(N)/b(N);
for j = (N-1):-1:1
y(j) = r(j)/b(j)-a(j)*y(j+1)/b(j);
end
plot(x,y)
4.8 二阶微分方程数值解
精确解为 $y(t) = \frac{1}{3}(t - \frac{1}{\sqrt{3}}\sin\sqrt{3}t)$。数值解的递推公式为 $y_{n + 1} = 2y_n - y_{n - 1} + \Delta t^2(-3y_n + t_n)$。
N = 20;
t = linspace(0,1,N);
dt = t(2)-t(1);
y(1) = 0;
y(2) = 0;
for j = 2:N-1
y(j+1) = 2*y(j)-y(j-1) + dt^2*(3*y(j)+t(j));
end
ex = (t-sin(sqrt(3)*t)/sqrt(3))/3;
plot(t,y,t,ex,'*')
4.9 二阶微分方程数值解
N = 20;
t = linspace(0,2,N);
dt = t(2)-t(1);
y(1) = 1;
y(2) = 1;
for j = 2:N-1
y(j+1) = 2*y(j)-y(j-1) + dt^2*(t(j)*y(j)+sin(t(j)));
end
plot(t,y)
另一个问题通过离散化方程 $y_{n + 1} - 2y_n + y_{n - 1} + \Delta t^2t_ny_n = \Delta t^2\sin t_n$ 并结合边界条件 $y_1 = 0$，$y_N = 0$ 求解。
N = 20;
t = linspace(0,2,N);
dt = t(2)-t(1);
A = zeros(N);
A = diag(-2/dt^2*ones(N,1)+t',0) + diag(1/dt^2*ones(N-1,1),-1) + diag(1/dt^2*ones(N-1,1),1);
r = transpose(sin(t));
A(1,:) = 0;
A(1,1) = 1; r(1) = 0;
A(N,:) = 0;
A(N,N) = 1; r(N) = 0;
sol = A\r;
plot(t,sol)
4.10 一阶微分方程
精确解为 $y(t) = \frac{1}{1 + t^2}$。数值解的递推公式为 $y_{n + 1} = y_n - \Delta t\frac{2t_ny_n}{1 + t_n^2}$。
N = 50;
t = linspace(0,5,N);
dt = t(2)-t(1);
y(1) = 1;
for j = 1:(N-1)
y(j+1) = y(j)-dt*2*t(j)*y(j)/(1+t(j)^2);
end
ex = 1./(1+t.^2);
plot(t,y,t,ex,'*')
4.11 线性系统的矩阵指数解
系统可表示为 $y’ = \begin{pmatrix}
0 & 1 & 0 \
0 & 0 & 1 \
2 & 1 & -2
\end{pmatrix}y$。通过计算矩阵的特征值和特征向量，得到 $e^{At} = Ve^{Dt}V^{-1}$，进而得到系统的解 $y(t) = e^{At}\begin{pmatrix}
0 \
0 \
1
\end{pmatrix}$。
A = [0 1 0; 0 0 1; 2 1 -2];
[V,D] = eig(A);
% 非归一化特征向量
V_tilde = [-1 1 1; 1 -2 1; -1 4 1];
V_tilde_inv = [-1 1/2 1/2; 1/2 -1/3 0; 1/3 1/3 1/6];
e_At = V_tilde * diag([exp(-t) exp(-2*t) exp(t)]) * V_tilde_inv;
y = e_At * [0; 0; 1];
disp('系统的解为 ');
disp(y);
4.12 有限差分法求解特征值问题
N = 20;
x = linspace(0,pi,N);
h = x(2)-x(1);
% Only the internal points
A = diag(-2/h^2*ones(N-2,1) + sin(x(2:(N-1)))',0) + diag(1/h^2*ones(N-3,1),1) + diag(1/h^2*ones(N-3,1),-1);
[V,D] = eigs(A,3,'SM');
for j = 1:3
ti = ['\lambda = ' num2str(D(j,j))];
subplot(1,3,j)
plot(x,[0; V(:,j); 0])
title(ti,'FontSize',14)
end
总结
本文涵盖了多项式拟合、插值、线性代数运算、积分计算和微分方程求解等多个方面的 MATLAB 实现。通过具体的代码示例和详细的解释，展示了如何使用 MATLAB 解决各种数值计算和线性代数问题。在实际应用中，可根据具体问题选择合适的方法和代码进行求解。
表格总结
| 问题类型 | 关键函数或方法 | 代码示例 |
|---|---|---|
| 多项式拟合 |
polyfit
|
p = polyfit(x, y, 1);
|
| 样条插值 |
spline
|
f = spline(x,y,z);
|
| 最小二乘法 | 矩阵运算 |
A = [sum(sin(x).^2) sum(cos(x).*sin(x)); sum(cos(x).*sin(x)) sum(cos(x).^2)]; r = [sum(f.*sin(x)); sum(f.*cos(x))]; sol = A\r;
|
| 线性方程组求解 |
\
|
x = A\b;
|
| 积分计算 |
sum
|
integral = (sum(f)-f(1)/2-f(N)/2)*h;
|
| 微分方程求解 |
ode45
|
[t,y] = ode45('odes',tspan,y0);
|
mermaid 流程图 - 多项式拟合流程
graph TD;
A[输入点坐标] --> B[构建向量 x 和 y];
B --> C[使用 polyfit 拟合直线];
C --> D[获取拟合参数 m 和 c];
D --> E[显示拟合结果];
mermaid 流程图 - 线性方程组求解流程
graph TD;
A[定义矩阵 A 和向量 b] --> B[使用 \ 运算符求解];
B --> C[得到解向量 x];
C --> D[显示解];
5. 综合应用与技巧总结
5.1 代码复用与模块化
在解决多个相关问题时，将常用的功能封装成函数可以提高代码的复用性和可维护性。例如，在积分计算中，将自定义函数封装成独立的
fn
函数：
function [val] = fn(x)
val = log(x+sqrt(x.^2+1));
end
这样在不同的积分计算场景中都可以方便地调用该函数。
5.2 误差分析与优化
在数值计算中，误差分析是非常重要的。例如，在插值和积分计算中，通过比较数值解和精确解可以评估误差。对于误差较大的情况，可以考虑增加采样点数或使用更精确的算法。
-
插值误差
：在插值误差计算中，通过比较
polyval(p,xv)
和
polyval(co,xv)
得到误差
err
，可以根据误差大小调整插值方法或增加数据点。
-
积分误差
：在不同点数的积分误差计算中，通过绘制误差与点数的关系图，观察误差的变化趋势，选择合适的点数进行积分计算。
5.3 矩阵运算的优化
在进行矩阵运算时，要注意矩阵的大小和运算的复杂度。例如，在矩阵乘法中，尽量避免不必要的重复计算。同时，利用 MATLAB 的内置函数可以提高运算效率。
-
矩阵乘法
：使用
A*B
进行矩阵乘法时，确保矩阵的维度匹配，避免出现维度不兼容的错误。
-
矩阵求逆
：在求解线性方程组
x = A\b
时，MATLAB 会根据矩阵的性质选择合适的算法，避免直接求逆矩阵，提高计算效率。
5.4 图形可视化
图形可视化可以帮助我们更直观地理解数据和结果。在 MATLAB 中，可以使用
plot
函数绘制各种图形。例如，在样条插值和积分误差分析中，通过绘制图形可以清晰地看到结果的变化趋势。
-
样条插值
：使用
plot(z,f,x,y,'o','MarkerSize',14)
绘制样条插值的结果，直观地展示插值曲线和原始数据点。
-
误差分析
：使用
plot(Ns,abs(integral-2))
绘制积分误差与点数的关系图，帮助我们选择合适的点数进行积分计算。
表格 - 代码优化技巧总结
| 优化类型 | 具体方法 | 示例代码 |
|---|---|---|
| 代码复用 | 封装函数 |
function [val] = fn(x); val = log(x+sqrt(x.^2+1)); end
|
| 误差分析 | 比较数值解和精确解 |
err = sum((true-f).^2);
|
| 矩阵运算优化 | 避免重复计算，利用内置函数 |
x = A\b;
|
| 图形可视化 |
使用
plot
函数
|
plot(z,f,x,y,'o','MarkerSize',14);
|
mermaid 流程图 - 代码优化流程
graph TD;
A[编写初始代码] --> B[进行误差分析];
B --> C{误差是否可接受};
C -- 是 --> D[结束优化];
C -- 否 --> E[优化代码];
E --> F[再次进行误差分析];
F --> C;
6. 实际案例分析
6.1 物理模型模拟
假设我们要模拟一个简单的物理系统，例如弹簧振子的运动。弹簧振子的运动可以用二阶微分方程来描述：
$m\frac{d^2x}{dt^2} + kx = 0$
其中 $m$ 是质量，$k$ 是弹簧的劲度系数。我们可以将其转化为一阶线性方程组：
$\begin{cases}
\frac{dx}{dt} = v \
\frac{dv}{dt} = -\frac{k}{m}x
\end{cases}$
使用 MATLAB 求解该方程组：
% 定义参数
m = 1;
k = 1;
% 定义微分方程
function dydt = spring_oscillator(t, y)
dydt = [y(2); -k/m * y(1)];
end
% 初始条件
y0 = [1; 0];
tspan = [0 10];
% 求解微分方程
[t, y] = ode45(@spring_oscillator, tspan, y0);
% 绘制结果
plot(t, y(:, 1), 'b-', t, y(:, 2), 'r--');
xlabel('时间');
ylabel('位置和速度');
legend('位置', '速度');
6.2 数据分析与预测
假设我们有一组数据，要对其进行分析和预测。例如，我们有一组股票价格数据，要预测未来的价格走势。我们可以使用多项式拟合或样条插值的方法对数据进行拟合，然后根据拟合结果进行预测。
% 生成示例数据
t = 1:10;
prices = [10, 12, 15, 13, 16, 18, 20, 19, 22, 25];
% 多项式拟合
p = polyfit(t, prices, 3);
% 预测未来价格
future_t = 11:15;
future_prices = polyval(p, future_t);
% 绘制结果
plot(t, prices, 'bo', future_t, future_prices, 'r*-');
xlabel('时间');
ylabel('股票价格');
legend('历史数据', '预测数据');
表格 - 实际案例总结
| 案例类型 | 问题描述 | 解决方法 | 代码示例 |
|---|---|---|---|
| 物理模型模拟 | 弹簧振子运动 |
转化为一阶线性方程组，使用
ode45
求解
|
[t, y] = ode45(@spring_oscillator, tspan, y0);
|
| 数据分析与预测 | 股票价格预测 | 多项式拟合，根据拟合结果预测 |
p = polyfit(t, prices, 3); future_prices = polyval(p, future_t);
|
mermaid 流程图 - 实际案例解决流程
graph TD;
A[问题描述] --> B[建立数学模型];
B --> C[选择合适的方法];
C --> D[编写代码求解];
D --> E[分析结果];
E --> F{结果是否满足要求};
F -- 是 --> G[结束];
F -- 否 --> B;
7. 学习建议与拓展
7.1 学习建议
- 多实践 ：通过实际编写代码来加深对 MATLAB 函数和方法的理解。可以从简单的问题开始，逐步增加难度。
- 阅读文档 ：MATLAB 官方文档是学习的重要资源，详细阅读文档可以了解函数的使用方法和参数含义。
- 参考示例 ：参考 MATLAB 提供的示例代码和教程，学习优秀的编程风格和解决问题的思路。
7.2 拓展学习
- 高级算法 ：学习更高级的数值计算算法，如遗传算法、粒子群算法等，用于解决更复杂的优化问题。
- 机器学习 ：结合 MATLAB 的机器学习工具箱，学习机器学习算法，如神经网络、支持向量机等，用于数据分析和预测。
- 图像处理 ：学习 MATLAB 的图像处理工具箱，进行图像的处理和分析，如图像滤波、特征提取等。
表格 - 学习资源总结
| 学习类型 | 资源推荐 |
|---|---|
| 官方文档 | MATLAB 官方网站的文档和教程 |
| 书籍 | 《MATLAB 数值计算与应用》等 |
| 在线课程 | Coursera、edX 等平台上的 MATLAB 相关课程 |
mermaid 流程图 - 学习流程
graph TD;
A[学习基础知识] --> B[实践项目];
B --> C{遇到问题};
C -- 是 --> D[查阅文档和资料];
D --> B;
C -- 否 --> E[拓展学习];
E --> F[学习高级算法和应用];
F --> B;
总结与展望
总结
本文全面介绍了 MATLAB 在数值计算和线性代数问题求解中的应用，涵盖了多项式拟合、插值、线性代数运算、积分计算、微分方程求解等多个方面。通过具体的代码示例和详细的解释，展示了如何使用 MATLAB 解决各种实际问题。同时，还介绍了代码优化、误差分析、实际案例应用等内容，帮助读者更好地掌握 MATLAB 的使用技巧。
展望
随着科学技术的不断发展，MATLAB 在各个领域的应用将越来越广泛。未来，我们可以进一步探索 MATLAB 在大数据分析、人工智能、深度学习等领域的应用，利用 MATLAB 的强大功能解决更复杂的问题。同时，不断学习和掌握新的算法和技术，提高自己的编程能力和解决问题的能力。
表格 - 整体内容总结
| 主题 | 主要内容 |
|---|---|
| 多项式拟合与插值 | 直线拟合、二次多项式拟合、样条插值等 |
| 线性代数问题 | 矩阵运算、线性方程组求解、特征值计算等 |
| 积分问题 | 积分计算、辛普森规则、弧长计算等 |
| 微分方程问题 | 一阶和二阶微分方程的精确解和数值解 |
| 综合应用与技巧 | 代码复用、误差分析、图形可视化等 |
| 实际案例分析 | 物理模型模拟、数据分析与预测 |
| 学习建议与拓展 | 学习方法、高级算法和应用的拓展 |
mermaid 流程图 - 整体知识体系
graph TD;
A[MATLAB 基础] --> B[数值计算];
A --> C[线性代数];
B --> D[多项式拟合与插值];
B --> E[积分计算];
B --> F[微分方程求解];
C --> G[矩阵运算];
C --> H[线性方程组求解];
C --> I[特征值计算];
D --> J[综合应用与技巧];
E --> J;
F --> J;
G --> J;
H --> J;
I --> J;
J --> K[实际案例分析];
K --> L[学习与拓展];
扫一扫
642
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
