求解变分不等式算例.doc

求解变分不等式：例2：（可以是维数，在我们计算的过程中，可以取100,200,1000维），是一个对称矩阵，可随机生成。例 1：其中,.其解为（4/3,7/9,4/9,2/9）。(不变)用迭代序列编程求解：高维迭代clc;k=0;k_inner=1000;time0=cputime;n=4;v0=0*rand(n,1);p0=1;Q=eye(n);b=5*diag(Q);%盒子的上界a=-5*diag(Q)%盒子的下届mu=0.03;%F函数的输入如下F=zeros(n,1)for i=1:nF(i,1)=v0(i)+exp(v0(i);endbarp0=max(0,p0-mu*sum(v0.2);% barv0=zeros(4,1);for i=1:n if v0(i)-mu*(F(i)+barp0*v0(i)=b(i) barv0(i)=b(i); else barv0(i)=v0(i)-mu*(F(i)+barp0*v0(i); endend while norm(v0-barv0)1e-5 | norm(p0-barp0)1e-5 & k=k_innerp0=max(p0-mu*sum(barv0.2);F=zeros(n,1)for i=1:nF(i,1)=barv0(i)+exp(barv0(i);endfor i=1:4 if v0(i)-mu*(F(i)+barp0*barv0(i)=b(i) v0(i)=b(i); else v0(i)=v0(i)-mu*(F(i)+barp0*barv0(i); endendF=zeros(n,1)for i=1:nF(i,1)=v0(i)+exp(v0(i);endbarp0=max(0,p0-mu*sum(v0.2);% barv0=zeros(4,1);for i=1:4 if v0(i)-mu*(F(i)+barp0*v0(i)=b(i) barv0(i)=b(i); else barv0(i)=v0(i)-mu*(F(i)+barp0*v0(i); endendk=k+1;endfprintf(vk is:%dn,v0);fprintf(p is:%dn,p0);fprintf(time used is:%dn,cputime-time0);fprintf(k is:%dn,k);迭代程序有限维clc;k=0;k_inner=1000;time0=cputime;v0=0;0;0;0;p0=1;A=4,2,2,1;2,4,0,1;2,0,2,2;-1,-1,-2,0;b=5;5;5;5;%盒子的上界a=-5;-5;-5;-5;%盒子的下届mu=0.03;q=-8;-6;-4;3;F=A*v0+q;barp0=max(0,p0-mu*(v0(1)2+v0(2)2+v0(3)2+v0(4)2);% barv0=zeros(4,1);for i=1:4 if v0(i)-mu*(F(i)+barp0*v0(i)=b(i) barv0(i)=b(i); else barv0(i)=v0(i)-mu*(F(i)+barp0*v0(i); endend while norm(v0-barv0)1e-5 | norm(p0-barp0)1e-5 & k=k_innerp0=max(p0-mu*(barv0(1)2+barv0(2)2+barv0(3)2+barv0(4)2);F=A*barv0+q;for i=1:4 if v0(i)-mu*(F(i)+barp0*barv0(i)=b(i) v0(i)=b(i); else v0(i)=v0(i)-mu*(F(i)+barp0*barv0(i); endendF=A*v0+q;barp0=max(0,p0-mu*(v0(1)2+v0(2)2+v0(3)2+v0(4)2);% barv0=zeros(4,1);for i=1:4 if v0(i)-mu*(F(i)+barp0*v0(i)=b(i) barv0(i)=b(i); else barv0(i)=v0(i)-mu*(F(i)+barp0*v0(i); endendk=k+1;用微分方程命令求解：function xdot=variational(t,x)A=4,2,2,1;2,4,0,1;2,0,2,2;-1,-1,-2,0;b=5;5;5;5;%盒子的上界a=-5;-5;-5;-5;%盒子的下届mu=0.03;q=-8;-6;-4;3;F=A*x(1);x(2);x(3);x(4)+q;%相当于变分不等式中的F(x)函数for i=1:4 if x(i)-mu*(F(i)+x(5)*x(i)=b(i) phi(i)=b(i)-x(i); else phi(i)=-mu*(F(i)+x(5)*x(i); endendxdot=phi(1);phi(2);phi(3);phi(4); max(0,x(5)+mu*(x(1)2+x(2)2+x(3)2+x(4)2)-x(5);t_final=100;x0=0;0;0;0;1;t,x=ode45(variational,0,t_final,x0);plot(t,x),figure;plot5(x(:,1),x(:,2),x(:,3),x(:,4),x(:,5);axis(10 40 -20 20 -20 20 );% %线性方程组的表达式：函数Gfunction G=variationalfunction(x)x=0;0;0;0;1;A=4,2,2,1;2,4,0,1;2,0,2,2;-1,-1,-2,0;b=5;5;5;5;%盒子的上界a=-5;-5;-5;-5;%盒子的下届mu=0.03;q=-8;-6;-4;3;F=A*x(1);x(2);x(3);x(4)+q;%相当于变分不等式中的F(x)函数for i=1:4 if x(i)-mu*(F(i)+x(5)*x(i)=b(i) phi(i)=b(i)-x(i); else phi(i)=-mu*(F(i)+x(5)*x(i); endendG=phi(1);phi(2);phi(3);phi(4); max(0,x(5)+mu*(x(1)2+x(2)2+x(3)2+x(4)2)-x(5);%函数G的Jacobian阵的计算：function Jacobian=Jacobian(x)x=0;0;0;0;1;A=4,2,2,1;2,4,0,1;2,0,2,2;-1,-1,-2,0;b=5;5;5;5;%盒子的上界a=-5;-5;-5;-5;%盒子的下届mu=0.03;Jacobian=zeros(5,5);q=-8;-6;-4;3;F=A*x(1);x(2);x(3);x(4)+q;%相当于变分不等式中的F(x)函数W=zeros(4,4);for i=1:4if x(i)-mu*(F(i)+x(5)*x(i)=b(i) W(i,i)=0; else W(i,i)=1; endendV=(eye(4)-mu*(A+x(5)*eye(4)*W-eye(4);for i=1:4;j=1:4;Jacobian(i,j)=V(i,j);endif x(5)+mu*(x(1)2+x(2)2+x(3)2+x(4)2)=0Jacobian(5,5)=1;elseJacobian(5,5)=0;end%x0=0;0;0;0;1;k=0;time0=cputime;G=variationalfunction(x0); while norm(G)1e-9 Jacobian=Jacobian(x0); d=inv(Jacobian)*G; x0=x0+d; k=k+1; G=variationalfunction(x0);endfprintf(xk is:%dn,x);fprintf(time used is:%dn,cputime-time0);fprintf(k is:%dn,k);%function xdot=variationalnew(t,x)A=4,2,2,1;2,4,0,1;2,0,2,2;-1,-1,-2,0;b=5;5;5;5;%盒子的上界a=-5;-5;-5;-5;%盒子的下届mu=0.03;B=4,8,8,8;9,5,0,5;4,2,6,7;4,6,3,4;g=x(1);x(2);x(3);x(4)*B*x(1);x(2);x(3);x(4);gradientg=(B+B)*x(1);x(2);x(3);x(4);q=-8;-6;-4;3;F=A*x(1);x(2);x(3);x(4)+q