数学与应用数学091
梁启凡
090080
实验一
实验目的：
编程并测试算法A.1 Gaussian Elimination with Row Partial Pivoting 和A.2 Cholesky Factorization
实验过程
一、A.1 Gaussian Elimination with Row Partial Pivoting
根据书上附录提供的伪代码，用MATLAB写出程序。用随机方阵测试，将测试结果与MATLAB自带的lu()函数产生的结果对比。
程序代码：
function [L,U,P]=GaussElimination(A)
%判断A是否奇异
if det(A)==0
error('matrix A is singular!')
stop
end
n=length(A);%求方阵A的维数
L=zeros(n);
P=eye(n);
for i=1:1:n
[m,t]=max(abs(A(i:n,i)));%找到A下三角阵中每列的最大值
j=t+i-1;%返回该最大值的行坐标
if(i~=j)
%交换A的i,j列
temp=A(i,1:n);
A(i,1:n)=A(j,1:n);
A(j,1:n)=temp;
%交换P的i,j列
temp=P(i,1:n);
P(i,1:n)=P(j,1:n);
P(j,1:n)=temp;
%交换L的i,j列
temp=L(i,1:n);
L(i,1:n)=L(j,1:n);
L(j,1:n)=temp;
end
L(i,i)=1;
for k=i+1:n
L(k,i)=A(k,i)/A(i,i);
for l=i+1:n
A(k,l)=A(k,l)-L(k,i)*A(i,l);
end
end
end
U=triu(A);%U为变化后的A的上三角部分
测试：
生成5维随机矩阵A：
>> A=rand(5)
A =
0.8147 0.0975 0.1576 0.1419 0.6557
0.9058 0.2785 0.9706 0.4218 0.0357
0.1270 0.5469 0.9572 0.9157 0.8491
0.9134 0.9575 0.4854 0.7922 0.9340
0.6324 0.99 0.8003 0.9595 0.6787
对A进行分解：
>> [L,U,P]=GaussElimination(A)
L =
1.0000 0 0 0 0
0.20 1.0000 0 0 0
0.1390 -0.5469 1.0000 0 0
0.9917 0.8870 0.9924 1.0000 0
0.6923 -0.3991 0.4794 0.1476 1.0000
U =
0.9134 0.9575 0.4854 0.7922 0.9340
0 -0.7565 -0.2753 -0.58 -0.1774
0 0 0.7391 0.4967 0.6223
0 0 0 -0.3559 -1.3507
0 0 0 0 -0.1376
P =
0 0 0 1 0
1 0 0 0 0
0 0 1 0 0
0 1 0 0 0
0 0 0 0 1
与MATLAB自带的lu函数结果完全相同：
二、A.2 Cholesky Factorization
根据书上附录提供的伪代码，用MATLAB写出程序。用随机方阵测试，将测试结果与MATLAB自带的chol()函数产生的结果对比。
程序代码：
function [L]=cholesky(A)
n=length(A);
for i=1:n
L(i,i)=sqrt(A(i,i));
for j=i+1:n
L(j,i)=A(j,i)/L(i,i);
for k=i+1:j
A(j,k)=A(j,k)-L(j,i)*L(k,i);
end
end
end
测试：
由于Cholesky分解要求被分解方阵是对称正定的，这里选取5阶帕斯卡矩阵进行测试。
即A。
用MATLAB自带函数chol()
与自己编写的分解函数所得结果互为转置。
实验分析：
以上两个程序均得到了正确结果，变换不同的维数的矩阵，也皆能得到正确的结果，说明程序是正确的。
对第二个程序，由于没有对称正定检测，所以哪怕输如不满足条件的矩阵也能得出结果。所以应当在输入A时确保A是对称正定的。
无约束最优化测试问题：
, ,
,
,
实验二
实验目的：
使用最速下降法，编程并求解。比较不同的终止参数和步长的选取。
实验过程及方案：
问题一的解决
终止参数选择：梯度的范数小于目标精确值
步长选择：固定步长 ，Armijo搜索
A：代码（固定步长）
function SDP1
syms x1 x2;
f=100*(x2-x1^2)^2+(1-x1)^2;
x0=[-1.2 1]';
eps=10e-3;
d1=jacobian(f);
d1=d1';
d2=jacobian(d1);
a=subs(d1,{x1,x2},{x0(1),x0(2)});
n=0;
while(norm(a)>eps)
p=-a;
step=0.0005;
x0=x0+step*p
a=subs(d1,{x1,x2},{x0(1),x0(2)});
n=n+1;
end
x0
n %迭代步数
结果：
x0 =
0.98
0.9779
n =
18179
迭代步长为0.0005，目标精度值为，得，共迭代了18179步，非常慢。
B:代码（Armijo搜索）
定义两个函数
%目标函数
function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1);
%目标函数的梯度，用jacobian（）函数求得
function gf=gfun(x)
gf=[ 2*x(1) - 400*x(1)*(- x(1)^2 + x(2)) - 2, - 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40, 2*x(3) - 360*x(3)*(- x(3)^2 + x(4)) - 2, - 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40]';
function SDP1(fun,gfun)
syms x
x0=[-1.2 1]';
maxk=5000000;%最大迭代步数
rho=0.5;
sigma=0.4;
k=0;
eps=10e-3;
while(k d=-g; if(norm(d) while(m<20)%Armoji搜索 if(feval(fun,x0+rho^m*d) end m=m+1; end x0=x0+rho^mk*d; k=k+1; end x0 %最优点 k %迭代次数 结果： x0 = 0.95 0.9791 k = 1202 达到同样的精度速度是很快的。 将精度改为10e-5 x0 = 1.0000 1.0000 k = 1435 精度相差10e2，搜索次数仅仅增加了200余次，效果拔群。 问题二的解决 终止参数选择：梯度的范数小于目标精确值 步长选择：Armijo搜索，wolfe线搜索条件 A：代码（Armijo搜索） 定义两个函数 %目标函数 function f=fun(x) f=100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1); %目标函数的梯度，用jacobian（）函数求得 function gf=gfun(x) gf=[ 2*x(1) - 400*x(1)*(- x(1)^2 + x(2)) - 2, - 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40, 2*x(3) - 360*x(3)*(- x(3)^2 + x(4)) - 2, - 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40]'; function SDP2(fun,gfun) syms x x0=[-3 -1 -3 -1]'; maxk=5000000;%最大迭代步数 rho=0.5; sigma=0.4; k=0; eps=10e-5; while(k d=-g; if(norm(d) while(m<20)%Armoji搜索 if(feval(fun,x0+rho^m*d) end m=m+1; end x0=x0+rho^mk*d; k=k+1; end x0 %最优点 k %迭代次数 结果： x0 = 1.0000 1.0001 1.0000 0.9999 k = 7156 非常精确。 B:代码（wolfe线搜索条件）fun函数和gfun函数不变 function SDP2(fun,gfun) syms x x0=[-3 -1 -3 -1]'; maxk=5000000;%最大迭代步数 c1=10e-4; c2=0.9; k=0; eps=10e-5; while(k d=-g; if(norm(d) step=1; while(step>0.0005)%wolfe搜索÷ if(feval(fun,x0+step*d) step=step;break; end step=step*0.8; end x0=x0+step*d k=k+1 end x0 %最优点 k %迭代次数 结果： x0 = 1.0000 1.0001 1.0000 0.9999 k = 21471 相同精度要求下，没有Armoji搜索快。 问题三的解决 终止参数选择：梯度的范数小于目标精确值 步长选择：Armijo搜索，wolfe线搜索条件 A：代码（Armijo搜索） 定义两个函数 %目标函数 function f=fun(x) f=(x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4; %目标函数的梯度，用jacobian（）函数求得 function gf=gfun(x) gf=[ 2*x(1) + 20*x(2) + 40*(x(1) - x(4))^3, 20*x(1) + 200*x(2) + 4*(x(2) - 2*x(3))^3, 10*x(3) - 10*x(4) - 8*(x(2) - 2*x(3))^3, 10*x(4) - 10*x(3) - 40*(x(1) - x(4))^3]’; function SDP3(fun,gfun) syms x x0=[3 -1 0 1]'; maxk=5000000;%最大迭代步数 rho=0.5; sigma=0.4; k=0; eps=10e-5; while(k d=-g; if(norm(d) while(m<20)%Armoji搜索 if(feval(fun,x0+rho^m*d) end m=m+1; end x0=x0+rho^mk*d; k=k+1; end x0 %最优点 k %迭代次数 结果： x0 = 0.0237 -0.0024 0.0118 0.0118 k = 13624 将精度eps改为10e-6 x0 = 0.0110 -0.0011 0.0055 0.0055 k = 62671 越来越接近实际值，但是搜索速度非常慢。 B:代码（wolfe线搜索条件）fun函数和gfun函数不变 代码同问题二的B方案一致。将代码中初始迭代点x0改为x0=[3 -1 0 1]' 结果： 速度非常慢，在进行到105416步的时候停止了计算，得到 x0 = 0.0468 -0.0047 0.0233 0.0234 这个结果离精确值差距还很大。 实验分析： 问题一用了2种步长，如果固定步长会使得搜索过程很慢，而逐一选择步长可以加快搜索过程和精度。问题二和问题三分别用Armoji搜索和wolfe搜索算步长的方法，Armoji搜索的速度要比wolfe快很多 实验三 实验目的： 使用牛顿法，编程并求解。 实验过程及方案： 问题一的解决 代码： function NTP1 syms x1 x2; f=100*(x2-x1^2)^2+(1-x1)^2; x0=[-1.2 1]'; eps=10e-3; d1=jacobian(f); d2=jacobian(d1); d1=d1'; a=subs(d1,{x1,x2},{x0(1),x0(2)}); b=subs(d2,{x1,x2},{x0(1),x0(2)}); n=0; while(norm(a)>eps) x0=x0-inv(b)*a a=subs(d1,{x1,x2},{x0(1),x0(2)}); b=subs(d2,{x1,x2},{x0(1),x0(2)}); n=n+1; end x0 n 结果： x0 = 1.0000 1.0000 n = 5 非常快速，只用5步就得到精度为10e-3的解。 问题二的解决 代码： function nt syms x1 x2 x3 x4; f=100*(x1^2-x2)^2+(x1-1)^2+(x3-1)^2+90*(x3^2-x4)^2+10.1*((x2-1)^2+(x4-1)^2)+19.8*(x2-1)*(x4-1); x0=[-3 -1 -3 -1]'; eps=10e-5; d1=jacobian(f); d2=jacobian(d1); d1=d1'; a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); n=0; while(norm(a)>eps) x0=x0-inv(b)*a a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); n=n+1; end x0 n 结果： x0 = -0.9680 0.9471 -0.9695 0.9512 n = 13 得到精度为10e-5的解只用了13步。 问题三的解决 代码： function nt syms x1 x2 x3 x4; f=(x1+10*x2)^2+5*(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4;; x0=[3 -1 0 1]'; eps=10e-20; d1=jacobian(f); d2=jacobian(d1); d1=d1'; a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); n=0; while(norm(a)>eps) x0=x0-inv(b)*a a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)}); n=n+1; end x0 n 结果： x0 = 1.0e-006 * 0.1435 -0.0144 0.0229 0.0229 n = 41 经过41步，得到精确度为10e-20的解。 实验分析： 牛顿法的速度较最速下降法有非常明显的改善，速度很快 实验四 实验目的： 调用MATLAB自带函数求解。 实验过程及方案： 分别用fminsearch函数和fminunc函数测试，对比两个函数的结果。 P1 1、使用fminsearch函数 >> f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x'); >> x0=[-1.2; 1]; ff=optimset; ff.Display='iter'; >> x=fminsearch(f,x0,ff) Iteration Func-count min f(x) Procedure 0 1 24.2 1 3 20.05 initial simplex 2 5 5.1618 expand 3 7 4.4978 reflect 4 9 4.4978 contract outside ……………………………………………………………………………… 81 151 1.12293e-009 contract inside 82 153 1.12293e-009 contract outside 83 155 1.12293e-009 contract inside 84 157 1.10755e-009 contract outside 85 159 8.17766e-010 contract inside Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 x = 1.0000 1.0000 总共用了85步得到结果。 2、使用fminunc函数 在刚才基础上输入 >>x=fminunc(f,[-1.2;.1],ff) First-order Iteration Func-count f(x) Step-size optimality 0 3 24.2 216 1 9 4.28049 0.000861873 15.2 2 12 4.12869 1 3 3 15 4.12069 1 1.5 …………………………………………………………………………………… 31 123 0.001075 0.5 0.877 32 126 9.14432e-005 1 0.143 33 129 7.27182e-006 1 0.0807 34 132 4.44847e-007 1 0.022 35 135 1.47502e-008 1 0.00272 36 138 2.83357e-011 1 1.91e-005 Local minimum found. Optimization completed because the size of the gradient is less than the default value of the function tolerance. x = 1.0000 1.0000 只用了36步便得到正确结果，比fminsearch函数快很多。 P2 1、使用fminsearch函数 >>f=inline('100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1)','x'); >> x0=[-3;-1;-3;-1]; ff=optimset; ff.Display='iter'; >> x=fminsearch(f,x0,ff) 314 527 1.94483e-009 contract inside 用了314步得到结果 x = 1.0000 1.0000 1.0000 1.0000 2、使用fminunc函数 >>x=fminunc(f,[-3;-1;-3;-1],ff) 35 215 1.84511e-007 1 0.00339 用35步得到结果 x = 1.0002 1.0004 0.9998 0.9996 结果没有fminsearch函数精确。 P3 1、使用fminsearch函数 >>f=inline('(x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4','x'); >> x0=[3;-1;0;1]; ff=optimset; ff.Display='iter'; >> x=fminsearch(f,x0,ff) 185 305 1.39059e-006 reflect 用185步得到结果 x = 0.0094 -0.0009 0.0166 0.0166 与真实值还有一些差距。 norm(x-[0 0 0 0]') =0.0253 2、使用fminunc函数 >>x=fminunc(f,[3;-1;0; 1],ff) 33 170 6.304e-009 1 0.000217 用33步得到结果 x = 0.0015 -0.0002 -0.0031 -0.0031 norm(x-[0 0 0 0]')=0.0047 比fminsearch函数得到的结果与精确值的差的范数要小。 实验分析： fminunc函数比fminisearch函数速度更快，快很多，不过精确性上是不分伯仲的。fminisearch与fminunc这两个函数都是最优化工具箱中的，如果需要速度快些，就用fminiunc函数。
Copyright © 2019-2026 51dongshi.net 版权所有
违法及侵权请联系：TEL:177 7030 7066 E-MAIL:11247931@qq.com 本站由北京市万商天勤律师事务所王兴未律师提供法律服务
