工程数学测试题（梯度下降、Gauss-Seidel迭代、拉格朗日插值）

1.手写程序，采用梯度下降算法，求最优解

解：函数接口，[x,result] = grad(f,x0)，其中f为目标函数，x0为初始值，x为最小值对应的坐标，result为迭代得到的最小值。

function [x,result]=grad(f,x0)
tol=1e-6;
f_sym=sym(f);
h=length(x0); 
FF='F_grad(';
ff='f(';
f_x1='f(';
for k=1:h
    FF=[FF,'x0(',num2str(k),'),'];
    ff=[ff,'x0(',num2str(k),'),'];
    f_x1=[f_x1,'x1(',num2str(k),'),'];
end
FF(end)=[];
FF=[FF,')'];
ff(end)=[];
ff=[ff,')'];
f_x1(end)=[];
f_x1=[f_x1,')'];
F_grad=matlabFunction(gradient(f_sym));
f_grad=eval(FF);
d_k=-f_grad/norm(f_grad);
while norm(f_grad) > tol
    syms alpha0
    x1=x0(:)+alpha0*d_k;
    fx1=eval(f_x1);
    d_x1=diff(fx1);
    d_alpha0=double(solve(d_x1));
    x0=x0(:)+d_alpha0*d_k;
    f_grad=eval(FF);
    if norm(f_grad) < tol
        break;
    end
    d_k=-f_grad/norm(f_grad);
end
x=x0;
result=eval(ff);

测试程序：

syms x1 x2
f=@(x1,x2)x1^2+6*x1*x2+10*x2^2+10*x1-4*x2;
[x,result]=grad(f,[0,0])

测试结果：

x =

  -56.0000
   17.0000

result =

 -314.0000

2.编写函数，采用Gauss-Seidel迭代法求解方程Ax=b

解：函数接口，result = gauss_seidel(A,b,x0)，其中，A 为系数矩阵，b 为常数向量，x0 为初始值；
输出result 为求解结果。

function result=gauss_seidel(A,b,x0)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)^(-1)*U;
f=(D-L)^(-1)*b;
flag=0;
k=0;
while(~flag)
    k=k+1;
   x=B*x0+f;
   if norm(x-x0)<1e-8
        result=x;
        flag=1;
   end
   if k>1000
       result=x;
       flag=1;
   end
   x0=x;
end

测试程序：

A=100*rand(3,3)
b=100*rand(3,1)
x0=[0;0;0];
result=gauss_seidel(A,b,x0)

测试结果：

A =

   91.7194   75.3729    7.5854
   28.5839   38.0446    5.3950
   75.7200   56.7822   53.0798


b =

   77.9167
   93.4011
   12.9906


result =

   -3.0999
    4.8595
   -0.5315

3.编写函数，进行拉格朗日插值

解：函数接口：p = lagrangeInterp(x, y)，其中，x是输入自变量序列，y是与每个自变量相应的函数值，输出p是多项式的系数序列。

function p=lagrangeInterp(x,y)
m=length(x);
for i=1:m
   p(i)=l(x,i)*y(i);
end
p=sym2poly(expand(sum(p)));
function L=l(x,n)
syms s
xn=x(n);
x(n)=[];
B=s-x;
C=xn-x;
L=prod(B)/prod(C);

测试程序：

x1=[2 3 5 7 8 ];
y1=[50 9 12 16 23 ];
p=lagrangeInterp(x1,y1)
x=x1(1):0.1:x1(end);
y=polyval(p,x);
plot(x,y)
xlabel('x')
ylabel('y')
title('LagrangeInterp')
hold on
plot(x1,y1,'o')

测试结果：

p =

0.5194  -11.6389   94.7139 -327.1944  410.3333