前言：
数值计算方法实验可以使用Matlab、C/C++、Python、Java等语言进行编程，考虑到同学期数学实验课程使用Matlab进行，建议提前熟悉Matlab编程（也效率更高）。
本文中各实验使用Matlab编程实现。
据观察，不同学年实验题目会有所变化，故本文将会介绍Matlab使用的一些基础知识，帮助理解。
（本文代码已添加注释）
目录
Matlab使用入门
实验一：非线性方程求根
共含六个脚本(.m)文件，实现二分法、简单迭代法、牛顿迭代法
f.m
%原函数为f=exp(x)+3*x^3-2*x^2+log(x)-1;
function f=fun(x)
f=exp(x)+3*x^3-2*x^2+log(x)-1;
end
g.m
%牛顿迭代法 x1=x-f(x)/f`(x) 公式对应函数
%原函数的导数为f'(x)=exp(x)+9*x^2-4*x+1.0/x
function g = g(x)
g=x-fun(x)/(exp(x)+9*x^2-4*x+1.0/x);
end
%或使用如下公式(可以直接求导数)
%syms a;
%f = a - (fun(a)./diff(fun(a)));
%g = subs(f,x);%牛顿迭代公式
main.m
%采用两种以上不同的算法求解exp(x)+3*x^3-2*x^2+ln(x)-1=0在[0,1]范围内的一个根，要求两次迭代误差小于1.0e-4。
%计算时将main.m中代码贴入命令行窗口
x0=0;
x1=0.5;
x2=1;
eps=1.0e-4;%计算精度
first_fun(x0,x2,eps);%二分法迭代
second_fun(x1,eps);%简单迭代法迭代
third_fun(x2,eps);%牛顿迭代法迭代
%plot:生成函数图像
x=linspace(0,1);
y=exp(x)+3*x.^3-2*x.^2+log(x)-1;
plot(x,y)
first_fun.m
%二分法计算函数的零点
%函数输出根的近似值和迭代次数
function [temp,k]=first_fun(down,up,eps)
k=0;
temp=0;
while abs(up-down)>eps %循环迭代，未达精度要求继续执行
if k>20
break
end
temp=(up+down)/2;
y3=fun(temp);
if y3==0
up=temp;
down=temp;
elseif fun(down)*y3<0
up=temp;
else
down=temp;
end
k=k+1;
%test
%disp('二分法迭代次数='),disp(k);
%disp('近似解=');double(temp)
end
disp('二分法迭代次数='),disp(k);
disp('近似解=');double(temp)
second_fun.m
%简单迭代法计算函数的零点
%函数输出根的近似值和迭代次数
function [k,x1]=second_fun(x0,eps)
k=1;
x1=fun(x0);
%test
%disp('简单迭代法迭代次数='),disp(k);
%disp('近似解=');double(x1)
while(abs(x1-x0)>eps)
if k>20
break
end
x0=x1;
x1=fun(x1);
k = k +1;
%test
%disp('简单迭代法迭代次数='),disp(k);
%disp('近似解=');double(x1)
end
disp('简单迭代法迭代次数='),disp(k);
disp('近似解=');
third_fun.m
%牛顿迭代法计算函数的零点
%函数输出根的近似值和迭代次数
function [k,x1]=third_fun(x0,eps)
k=1;
x1=g(x0);
%test
%disp('牛顿迭代法迭代次数='),disp(k);
%disp('近似解='),disp(x1);
while abs(x0-x1)>eps %eval(x0-x1)
x0=x1;
if k>20
break
end
x1 =g(x0);
k = k+1;
%test
%disp('牛顿迭代法迭代次数='),disp(k);
%disp('近似解='),disp(x1);
end
disp('牛顿迭代法迭代次数='),disp(k);
disp('近似解='),disp(x1);实验二：线性方程组求解
共含五个脚本(.m)文件，实现高斯消元法、列主元消元法、雅可比迭代法、高斯-赛德尔迭代法
main.m
%已知如下四元一次线性方程组,请选择至少两种算法编程求解它，并分析所得到结果的合理性。
a=[14,2,1,5;8,17,2,10;4,18,3,6;12,26,11,20];%原系数矩阵
b=[1;2;3;4];%常数列向量
x0=[0;0;0;0];
eps=1.0e-6;%精度
%disp('高斯消元法')
c=gauss(a,b);
%disp('列主元消元法')
d=lzy(a,b);
%disp('Jacobi迭代法')
[e,n1]=Jacobi(a,b,x0,eps);
%disp('高斯赛德尔迭代法')
[f,n2]=gauss_seidel(a,b,x0,eps);
gauss.m
%高斯消元法求解
function x=gauss(a,b)
%参量说明：A为系数矩阵；B为常数列向量；extend为增广矩阵
%将增广矩阵化为上三角，再回带求解x
%此方法较为常规，将extend(k,k)元素乘以-extend(i,k)/extend(k,k)加到第i行
%从1:n-1列，主对角元素的以下行，通过两层循环来遍历
extend=[a,b];
n=length(b);
ra=rank(a);
rz=rank(extend);
temp=rz-ra;
if temp>0
disp('无解');
return;
end
if ra==rz
if ra==n
x=zeros(n,1);
for p=1:n-1
for k=p+1:n
m=extend(k,p)/extend(p,p);
extend(k,p:n+1)=extend(k,p:n+1)-m*extend(p,p:n+1);
end
end
b=extend(1:n,n+1);
a=extend(1:n,1:n);
x(n)=b(n)/a(n,n);
for q=n-1:-1:1
x(q)=(b(q)-sum(a(q,q+1:n)*x(q+1:n)))/a(q,q);
end
end
end
disp('高斯消元法解得:'),disp(x);
lzy.m
%列主元消元法求解
function x=lzy(a,b)
matrix=[a,b];
[row,col]=size(matrix);
x=zeros(row,1);
%判定输入矩阵是否符合要求
if row~=col-1
disp('请输入n*(n+1)维矩阵');
else
for ii = 1:row-1
%寻找主元
max = ii;
for jj = ii+1:row
if abs(matrix(jj,ii)) > abs(matrix(max,ii))
max = jj;
end
end
matrix([ii,max],:) = matrix([max,ii],:);
if matrix(ii,ii) == 0
disp(['第',num2str(ii),'个主元素为零']);
return;
end
%开始消元
for jj = ii+1:row
matrix(jj,:) = matrix(jj,:)-...
matrix(jj,ii)/matrix(ii,ii)*matrix(ii,:);
end
end
%消元完毕，开始回代
if matrix(row,row)==0
disp(['第',num2str(row),'个主元素为零']);
return;
end
x(row)=matrix(row,col)/matrix(row,col-1);
for ii=row-1:-1:1
x(ii)=(matrix(ii,col)...
-matrix(ii,1:row)*x)/matrix(ii,ii);
end
end
disp('列主元消元法解得:'),disp(x);
end
Jacobi.m
%Jacobi迭代法
function [y,n]=Jacobi(a,b,x0,eps)
D=diag(diag(a)); %求对角矩阵
L=-tril(a,-1); %求下三角矩阵
U=-triu(a,1); %求上三角矩阵
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=eps %用2范数去逼近
x0=y;
y=B*x0+f;
n=n+1;
end
disp('Jacobi迭代法解得:'),disp(y);
disp('迭代次数：'),disp(n);
gauss_seidel.m
%高斯赛德尔迭代法
function [y,n]=gauss_seidel(a,b,x0,eps)
D=diag(diag(a)); %求对角矩阵
L=-tril(a,-1); %求下三角阵
U=-triu(a,1); %求上三角阵
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=eps %用2范数去逼近
x0=y;
y=B*x0+f;
n=n+1;
end
disp('GaussSeidel迭代法解得:'),disp(y);
disp('迭代次数：'),disp(n);
实验三：数值积分
共含四个脚本(.m)文件，实现复合梯形公式、复合Simpson公式
fun.m
%原函数
function f = fun(x)
f=(2*x-x^2+x^3+exp(x/2))^0.5;
end
main.m
%选择一种或多种数值积分方法求解下列函数在区间[0,2]内的积分值，要求精度小于1.0e-5：
%f(x) =(2*x-x^2+x^3+exp(x/2))^0.5
a=0;
b=2;
e=1.0e-5;%精度
n=1;
Simpson(a,b,n,e);
tixing(a,b,n,e) ;
%plot:生成函数图像
%x=linspace(0,2);
%y=(2*x-x.^2+x.^3+exp(x/2)).^0.5;
%plot(x,y)
tixing.m
function [h,n] = tixing(a,b,n,e) %复合梯形公式
%区间为[a,b],n为初始区间切割次数,e为精度
%先求初始步长对应的公式值
h = (b-a)/n; %步长
T1= fun(a)+fun(b);
T1=T1*h/2;
%求变步长对应的公式值
n=n+1;
h = (b-a)/n;
x=a;
T2=fun(x)+fun(b);
for k=1:n-1
x = x+h;
T2 = T2+ 2*fun(x);
end
T2=T2*h/2;
while abs(T2-T1)>e %循环迭代，未达精度要求继续执行
T1=T2;
x=a;
T2=fun(x)+fun(b);
n=n+1;
h = (b-a)/n;
for k=1:n-1
x = x+h;
T2 = T2+ 2*fun(x);
end
T2=T2*h/2;
end
disp('复合梯形公式迭代积分值:'),disp(T2);
disp('最终步长:'),disp(h);
disp('迭代次数:'),disp(n);
Simpson.m
function [h,n] = Simpson(a,b,n,e) %复合Simpon公式
%区间为[a,b],n为初始区间切割次数,e为精度
%先求初始步长对应的公式值
h = (b-a)/n; %步长
x = a;
S1= fun(x)-fun(b);
for k=1:n
x = x+h/2;
S1 = S1+ 4*fun(x);
x=x+h/2;
S1=S1+ 2*fun(x);
end
S1=S1*h/6;
%求变步长对应的公式值
n=n+1;
h = (b-a)/n;
x=a;
S2=fun(x)-fun(b);
for k=1:n
x = x+h/2;
S2 = S2+ 4*fun(x);
x=x+h/2;
S2 = S2+ 2*fun(x);
end
S2=S2*h/6;
while abs(S2-S1)>e %循环迭代，未达精度要求继续执行
S1=S2;
x=a;
S2=fun(x)-fun(b);
n=n+1;
h = (b-a)/n;
for k=1:n
x = x+h/2;
S2 = S2+ 4*fun(x);
x=x+h/2;
S2 = S2+ 2*fun(x);
end
S2=S2*h/6;
end
disp('复合Simpon公式迭代积分值:'),disp(S2);
disp('最终步长:'),disp(h);
disp('迭代次数:'),disp(n);
