共轭梯度法C语言源程序
共轭梯度法的C+=语言实现
共轭梯度法的C+=语言实现---------------------基于成都理工大学最优化教材#include <iostream.h>#include "Matrix.h"#include<LIMITS.H>#define MAX_M 2048#include<math.h>#define beta (sqrt(5)-1)/2#define eclipse 0.01CMatrix hesse(){CMatrix temp(2,2);temp.A[0][0] = 2;temp.A[0][1] = 0;temp.A[1][0] = 0;temp.A[1][1] = 8;return temp;}double fun(double t,CMatrix &X1,CMatrix &X2){double x1 = X1.A[0][0] + t*X2.A[0][0];double x2 = X1.A[1][0] + t*X2.A[1][0];return x1*x1 + 4*x2*x2;}double Max_Num(double a,double b){if(a>b){return a;}else{return b;}}double Min_Num(double a,double b){if(a<b){return a;}else{return b;}}void StepAdding_Search(double &a,double &b,CMatrix &mt1,CMatrix &mt2) {double t[MAX_M]={0};double h[MAX_M]={0};double f[MAX_M]={0};double result=0;double p=2.1;t[0]=0;h[0]=1;f[0]=fun(t[0],mt1,mt2);for(int k=0;k<MAX_M-1;k++){t[k+1]=t[k]+h[k];f[k+1]=fun(t[k+1],mt1,mt2);if(f[k+1]<f[k]){h[k+1]=p*h[k];result=t[k];t[k]=t[k+1];f[k]=f[k+1];}else if(k == 0){h[k]=-h[k];result=t[k+1];}else{a=Min_Num(result,t[k+1]);b=Max_Num(result,t[k+1]);}}}double Fibonacci(double &a,double &b,CMatrix mt1,CMatrix mt2){double t2,t1,result_1,result_2;t2 = a + beta*(b-a);result_2 = fun(t2,mt1,mt1);while(true){t1 = a+b-t2;result_1 = fun(t1,mt1,mt2);if(fabs(t1-t2)<eclipse){return (t1+t2)/2;}else if(result_1 < result_2){b = t1;t1 = t2;result_1 = result_2;t2 = a+beta*(b-a);result_2 = fun(t2,mt1,mt2);}else if(result_1 == result_2){a = t1;b = t2;result_2 = fun(a,mt1,mt2);}else{a = t2;t2 = t1;result_2 = result_1;}}}double distance(CMatrix &X1,CMatrix &X2) {double x1 = X1.A[0][0] - X2.A[0][0];double x2 = X1.A[1][0] - X2.A[1][0];return x1*x1 + x2*x2;}CMatrix diff_fun(CMatrix &mt){CMatrix temp(2,1);temp.A[0][0] = 2*mt.A[0][0];temp.A[1][0] = 8*mt.A[1][0];return temp;}double fanshu(CMatrix mt){return sqrt(mt.A[0][0]*mt.A[0][0] + mt.A[1][0]*mt.A[1][0]);}void main(){int i =0;int n = 2000;CMatrix temp_X1;CMatrix X1(2,1);X1.A[0][0] = 1;X1.A[1][0] = 1;CMatrix m_temp = hesse();CMatrix n_temp = diff_fun(X1);CMatrix X2 = (m_temp.Inverse()* n_temp).Matrix_Multiply(-1);double a,b;StepAdding_Search(a,b,X1,X2);double m = Fibonacci(a,b,X1,X2);temp_X1 = X1;CMatrix X3 = X1 + X2.Matrix_Multiply(m);while(distance(X1,X3) > eclipse){X1 = X3;if(i+1 == n){i = 0;X2 = m_temp.Inverse()* diff_fun(X1);}else{m = (fanshu(temp_X1)*fanshu(temp_X1))/(fanshu(X1)*fanshu(X1));X2 = diff_fun(X1).Matrix_Multiply(-1) + X2.Matrix_Multiply(m);}X1 = X3;X2 = m_temp.Inverse()* diff_fun(X1);StepAdding_Search(a,b,X1,X2);m = Fibonacci(a,b,X1,X2);X3 = X1 + X2.Matrix_Multiply(m);i++;}cout<<"求解的结果是:"<<endl;cout<<"("<<X1.A[0][0]<<","<<X1.A[1][0]<<")"<<endl;cout<<"最优解是:"<<endl;cout<<fun(m,X1,X2)<<endl;}（注：本资料素材和资料部分来自网络，仅供参考。
共轭梯度法步骤
共轭梯度法步骤共轭梯度法是一种求解线性方程组的迭代算法，它以高效稳定的特点而广受欢迎。
以下是共轭梯度法的步骤：步骤1：初始化首先，我们需要有一个初始向量x0和一个初始残量r0=b-Ax0。
其中，A为系数矩阵，b为常数向量。
步骤2：计算方向向量令d0=r0，表示第一次迭代的方向向量。
步骤3：计算步进长度令α0=(r0·r0)/(d0·Ad0)，其中·表示向量的点积。
α0表示迭代过程中每个方向向量的步进长度。
步骤4：更新解向量令x1=x0+α0d0，表示迭代后的解向量。
步骤5：计算新残量令r1=r0-α0Ad0。
步骤6：判断终止条件如果r1的范数小于预设阈值，或者迭代次数达到预设次数，终止迭代。
否则，进入下一次迭代。
步骤7：更新方向向量令β1=(r1·r1)/(r0·r0)，表示更新方向向量的轴线。
步骤8：计算新方向向量令d1=r1+β1d0，表示新的迭代方向向量。
步骤9：计算新的步进长度令α1=(r1·r1)/(d1·Ad1)。
步骤10：更新解向量令x2=x1+α1d1。
步骤11：更新残量令r2=r1-α1Ad1。
步骤12：重复步骤6至11，直至满足终止条件。
总结起来，共轭梯度法的步骤主要包括初始化、计算方向向量、计算步进长度、更新解向量、计算新残量、判断终止条件、更新方向向量、计算新的步进长度、更新解向量和更新残量等。
该算法迭代次数较少，收敛速度快，适用于大规模线性方程组的求解。
共轭梯度法C语言(西安交大)
#include<stdio.h>#include<math.h>#define N 10 /*定义矩阵阶数*/void main(){int i,j,m,A[N][N],B[N];double X[N],akv[N],dka[N],rk[N],dk[N],pk,pkk,ak,bk;for(i=0;i<N;i++) /*输入系数矩阵A*/ {for(j=0;j<N;j++){if(i==j)A[i][j]=-2;else if(abs(i-j)==1)A[i][j]=1;elseA[i][j]=0;}}for(i=0;i<N;i++) /*输入矩阵B*/ {if(i==0||i==N-1)B[i]=-1;elseB[i]=0;}printf("\n"); /*输出系数矩阵A*/printf("the Maxtr A is\n");for(i=0;i<N;i++){printf("\n");for(j=0;j<N;j++)printf("%3d",A[i][j]);}printf("\n"); /*输出矩阵B*/printf("the Maxtr b is\n");for(i=0;i<N;i++)printf("%3d",B[i]);printf("\n");printf("input the Maxtr X\n"); /*给X0输入一个初始向量*/ for(i=0;i<N;i++)X[i]=0;printf("X chushi:\n");for(i=0;i<N;i++) /*显示给定的初始向量X0*/ printf("%3.0f",X[i]);/*开始计算*/for(i=0;i<N;i++) /*计算rk*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*X[j];akv[i]=pk;}for(i=0;i<N;i++)rk[i]=B[i]-akv[i];for(i=0;i<N;i++) /*给rO赋值*/dk[i]=rk[i];for(m=0;m<N;m++) /*开始迭代循环*/{for(i=0;i<N;i++) /*计算ak*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*dk[j];dka[i]=pk;}pk=0.0;pkk=0.0;for(i=0;i<N;i++){pk+=dka[i]*dk[i];pkk+=rk[i]*rk[i];}if(pkk==0) /*如果残差pkk=0,计算结束*/ break;ak=pkk/pk;for(i=0;i<N;i++) /*计算X(k+1)*/X[i]=X[i]+ak*dk[i];for(i=0;i<N;i++) /*计算r(k+1)*/{pk=0.0;for(j=0;j<N;j++)pk+=A[i][j]*X[j];akv[i]=pk;}for(i=0;i<N;i++)rk[i]=B[i]-akv[i];pk=0.0;for(i=0;i<N;i++) /*计算bk*/pk+=rk[i]*rk[i];bk=pk/pkk;for(i=0;i<N;i++) /*计算d(k+1)*/dk[i]=rk[i]+bk*dk[i];}printf("\n");printf("X cacualtate value is:\n"); /*输出运算结果X*/for(i=0;i<N;i++)printf("%f\n",X[i]);}。
共轭梯度法解油藏数值模拟线性方程中的稀疏矩阵及C++程序(测试通过)
for(i=0;i<=n-1;i++) {
x[i]=0.0; p[i]=b[i]; r[i]=b[i]; } i=0; while(i<=n-1) { for(k=0;k<=n-1;k++) {
s[k]=0.0; for(j=0;j<=n-1;j++)
s[k]=s[k]+a[k*n+j]*p[j]; } d=0.0; e=0.0; for(k=0;k<=n-1;k++) {
Байду номын сангаасfor(int jj=0;jj<5;jj++) {
aa[m++]=a[ii][jj];
} } double b[5]={23,32,33,31,66}; eps=0.000001; grad(5,aa,b,eps,x); for(i=0;i<=4;i++) {
cout<<"x["<<i<<"]="<<x[i]<<endl; } }
d=d+p[k]*b[k]; e=e+p[k]*s[k]; } alpha=d/e; for(k=0;k<=n-1;k++) { x[k]=x[k]+alpha*p[k];
} for(k=0;k<=n-1;k++) {
q[k]=0.0; for(j=0;j<=n-1;j++) {
q[k]=q[k]+a[k*n+j]*x[j]; } } d=0.0; for(k=0;k<=n-1;k++) { r[k]=b[k]-q[k]; d=d+r[k]*s[k]; } beta=d/e; d=0.0; for(k=0;k<=n-1;k++) { d=d+r[k]*r[k]; } d=sqrt(d); if(d<eps) {
共轭梯度法的C实现
西安交通大学实验报告课程名称：数值分析上机实验实验名称：共轭梯度法学院：___数学学院______________________ 班级姓名：学号：实验日期 2015 年 05 月 26 日自评成绩：97一、实验目的(1)熟练掌握改进平方根法和共轭梯度法的迭代过程(2)尝试使用自己熟悉的计算机语言解决数学中的问题(3)通过上机实验来巩固课本中所学的知识二、实验内容与结果题目2：共轭梯度法源程序2#include <iostream>using namespace std;double f1(double a[10],double n)//构造第一个求和函数，简化主函数{double s=0;int i;for(i=0;i<n;i++){s=s+a[i]*a[i];}return s;}double f2(double r[10],double a[10][10],int n)//构造第二个求和函数，简化主函数{double b[10],m=0;int i,j;for(i=0;i<n;i++){double s=0;for(j=0;j<n;j++){s=s+r[j]*a[j][i];}b[i]=s;}for(i=0;i<n;i++){m=m+b[i]*r[i];}return m;}double *f3(double a[10][10],double b[10],int n)//构造输出列向量的函数{double *r;r=new double[10];int i,j;for(i=0;i<n;i++){double s=0;for(j=0;j<n;j++){s=s+a[i][j]*b[j];}r[i]=s;}return r;}int main(){int i,j,k,n;double a,b,e,A[10][10],B[10],r[10],r1[10],x[10],d[10];double *ax;cout<<"请输入未知数的个数和计算精度："<<endl; cin>>n>>e;cout<<"请输入起始向量"<<endl;for(i=0;i<n;i++){cin>>x[i];}cout<<"请输入右端项： "<<endl;for(i=0;i<n;i++){cin>>B[i];}cout<<"请输入矩阵： "<<endl;for(i=0;i<n;i++){for(j=0;j<n;j++){cin>>A[i][j];}}ax=f3(A,x,n);for(i=0;i<n;i++) {r[i]=B[i]-ax[i];d[i]=r[i];}if(sqrt(f1(r,n))<=e) {for(i=0;i<n;i++) {cout<<x[i]<<'\t';}}else{for(k=0;k<n;k++){a=f1(r,n)/f2(d,A,n);for(i=0;i<n;i++){x[i]=x[i]+a*d[i];}ax=f3(A,x,n);for(i=0;i<n;i++){r1[i]=B[i]-ax[i];}if(sqrt(f1(r1,n))<=e||k+1==n) {for(i=0;i<n;i++){cout<<"x["<<i<<"]="<<x[i]<<'\n'; }break;}else{b=f1(r1,n)/f1(r,n);for(i=0;i<n;i++){d[i]=r1[i]+b*d[i];r[i]=r1[i];}}}}return 0;}运行结果2三、实验总结（列出实验的结论、收获和存在的问题，必写）通过这次实验，我对简化平方根法和共轭梯度法的计算过程和迭代格式有了熟练地掌握，让我学习到了很多，同时也是一个复习的过程。
共轭梯度（ConjugateGradient，CG）算法
共轭梯度（ConjugateGradient，CG）算法
以下皆为从⽹络资料获取的感性认知
共轭定义
共轭在数学、物理、化学、地理等学科中都有出现。
本意：两头⽜背上的架⼦称为轭，轭使两头⽜同步⾏⾛。
共轭即为按⼀定的规律相配的⼀对。
通俗点说就是孪⽣。
在数学中有共轭复数、共轭根式、共轭双曲线、共轭矩阵等。
求解⽬标
共轭梯度法是数值分析⾥⾯⼀类⾮常重要的⽅法，是⽤来解决⽆约束凸⼆次规划问题的⼀种⽅法
⽅法的⽬的就是通过迭代求得多元⼆次函数 [公式] 的极值点
基本思想
每个⽅向上只找⼀次，对于 [公式] 维空间就进⾏ [公式] 次计算，也就是说，我每个⽅向上都能⾛到极致，这样我就不会⾛任何的回头路了，效率也就最⾼。
那么是如何实现的呢？如果每个⽅向上都是正交的，我肯定就是在这个⽅向上⾛到极致了，这样⽅向也就确定了，也就是每⼀次都⾛正交⽅向，然后再确定每次的步长就可以了。
抽象图。
共轭梯度法的C程序代码
(1-x[1])*(1-x[1])+(1-x[4])*(1-x[4])+(x[1]*x[1]-x[2])*(x[1]*x[1]x
[2])+(x[2]*x[2]-x[3])*(x[2]*x[2]-x[3])+(x[3]*x[3]-x[4])*(x[3]*x[ 3]-x[4]);//例二
} double vecFun_Si(vector vx){ //不等式约束函数 double x[SIZE]; for(int i=0;i<vx.dimension;++i) x[i+1]=vx.elem[i]; return x[1]+1;//不等式约束函数 } double vecFun_Hi(vector vx){
double vecMultiply(vector a, vector b){ //行向量与列向量乘法运算 if((a.tag!=b.tag)&&(a.dimension==b.dimension)){//相乘条件 double c=0; for(int i=0;i<a.dimension;++i) c+=a.elem[i]*b.elem[i]; return c; } } vector vecAdd(vector a, vector b){ //向量加法运算 if((a.tag==b.tag)&&(a.dimension==b.dimension)){//相加条件 vector c; c.dimension=a.dimension; c.tag=a.tag; for(int i=0;i<c.dimension;++i) c.elem[i]=a.elem[i]+b.elem[i]; return c; } } vector vecConvert(vector a){ //向量转置 if(a.tag==0)a.tag=1;
共轭梯度法 c++
共轭梯度法c++一、共轭梯度法是一种优化算法，特别适用于解决对称正定矩阵的线性方程组。
它通过迭代的方式逐步逼近方程组的解，具有较快的收敛速度。
在C++中实现共轭梯度法可以为解决大规模线性系统提供高效的数值解。
二、共轭梯度法基本原理问题背景：考虑一个线性方程组Ax = b，其中A是对称正定矩阵，b是已知向量。
迭代过程：共轭梯度法通过迭代寻找一个逼近解x_k，使得残差r_k = b - Ax_k 最小。
迭代过程中，每一步都保证搜索方向共轭于前一步的搜索方向。
算法步骤：初始化：选择初始解x_0，计算残差r_0 = b - Ax_0，初始化搜索方向p_0 = r_0。
迭代：对于每一步k，计算步长alpha_k，更新解x_k+1 = x_k + alpha_k * p_k，计算新的残差r_k+1，更新搜索方向p_k+1。
收敛检测：当残差足够小时，停止迭代。
共轭方向的选择：在每一步中，选择共轭搜索方向可以通过Gram-Schmidt正交化方法得到。
这样能够确保搜索方向之间是线性无关的。
三、C++中的共轭梯度法实现在C++中实现共轭梯度法需要考虑以下关键步骤：矩阵和向量表示：使用C++中的数组或矩阵库表示矩阵A和向量b。
迭代过程：实现共轭梯度法的迭代过程，包括更新解、计算残差、计算步长等。
共轭方向选择：使用Gram-Schmidt正交化方法确保搜索方向共轭。
收敛检测：制定合适的收敛准则，如残差的阈值，判断是否停止迭代。
以下是一个简化的C++示例代码，演示了共轭梯度法的基本实现：#include <iostream>#include <cmath>#include <vector>using namespace std;// 定义矩阵和向量类型typedef vector<vector<double>> Matrix;typedef vector<double> Vector;// 共轭梯度法实现Vector conjugateGradient(const Matrix& A, const Vector& b, const Vector& x0, double tolerance, int maxIterations) {int n = A.size();Vector x = x0;Vector r = b - multiply(A, x);Vector p = r;for (int k = 0; k < maxIterations; ++k) {double alpha = dot(r, r) / dot(p, multiply(A, p));x = x + alpha * p;Vector newR = r - alpha * multiply(A, p);double beta = dot(newR, newR) / dot(r, r);p = newR + beta * p;r = newR;// 收敛检测if (sqrt(dot(r, r)) < tolerance) {cout << "Converged after " << k + 1 << " iterations." << endl;break;}}return x;}// 辅助函数：向量点积double dot(const Vector& a, const Vector& b) {double result = 0.0;for (size_t i = 0; i < a.size(); ++i) {result += a[i] * b[i];}return result;}// 辅助函数：矩阵与向量相乘Vector multiply(const Matrix& A, const Vector& x) {int n = A.size();Vector result(n, 0.0);for (int i = 0; i < n; ++i) {for (int j = 0; j < n; ++j) {result[i] += A[i][j] * x[j];}}return result;}int main() {// 示例用例Matrix A = {{4, 1}, {1, 3}};Vector b = {1, 2};Vector x0 = {0, 0};double tolerance = 1e-6;int maxIterations = 1000;// 调用共轭梯度法Vector solution = conjugateGradient(A, b, x0, tolerance, maxIterations);// 输出结果cout << "Solution: ";for (double value : solution) {cout << value << " ";}cout << endl;return 0;}这个简单的示例演示了如何使用C++实现共轭梯度法。
最优化问题共轭梯度法法代码
最优化问题共轭梯度法法代码x本文介绍了最优化问题共轭梯度法法的代码实现，以及如何使用代码解决最优化问题。
一、共轭梯度法简介共轭梯度法是一种常用的最优化算法，它是一种经典的迭代方法，用于求解凸函数的极值问题。
其基本思想是：在每一步，沿着梯度下降的方向迭代，直到梯度为零就停止迭代。
共轭梯度法的迭代公式为：$$x_{k+1}=x_k+alpha_k p_k$$其中，$alpha_k$ 是步长参数，$p_k$ 是当前搜索方向，$x_k$ 是当前点的位置。
二、代码实现1.函数定义```python# 共轭梯度法# 入参：函数func，梯度函数grad，初始点x0，步长参数alpha，精度epsilon# 出参：求解的最优点xdef conjugate_gradient_method(func, grad, x0, alpha, epsilon):```2.初始化搜索方向```python# 初始化搜索方向p_k = -grad(x_k)```3.更新迭代点```python# 更新迭代点x_k = x_k + alpha * p_k```4.更新搜索方向```python# 更新搜索方向beta_k = (grad(x_k) * grad(x_k)) / (grad(x_k_prev) * grad(x_k_prev))p_k = -grad(x_k) + beta_k * p_k_prev```5.检查终止条件```python# 检查终止条件if np.linalg.norm(grad(x_k)) < epsilon:break```6.完整代码```python# 共轭梯度法# 入参：函数func，梯度函数grad，初始点x0，步长参数alpha，精度epsilon# 出参：求解的最优点xdef conjugate_gradient_method(func, grad, x0, alpha, epsilon):x_k = x0p_k = -grad(x_k)while np.linalg.norm(grad(x_k)) > epsilon:x_k_prev = x_kp_k_prev = p_kline_search = line_search_method(func, grad, p_k, x_k, alpha)x_k = x_k + line_search * p_kbeta_k = (grad(x_k) * grad(x_k)) / (grad(x_k_prev) * grad(x_k_prev))p_k = -grad(x_k) + beta_k * p_k_prevreturn x_k```三、如何使用代码解决最优化问题1.确定问题首先，我们需要确定最优化问题，即构造一个函数，其中包含我们想要优化的目标函数以及约束条件。
梯度法搜索优化c语言源程序
void Get_Search_area(int n,double *x_k,double *s_k,double *a1,double *a3)
{
}
}
}
/*黄金分割法求极小值*/
/*入口参数：
n :优化模型维数
a1 :搜索区间步长下限
a4 :索区间步长上限
s_k :n维数组，搜索方向向量 */
/*出口参数：
x_k :n维数组，极小点坐标 */
}
/* 计算X(k+1)=X(k)+a*S(k)
入口参数：
x_k : n维数组，自变量
s_k : n维数组，方向向量
a : 方向步长
n : 维数
出口参数：
x_k1 : 转换后的n维数组*/
void fun(double *x_k1,double *x_k,double *s_k,double a,int n)
a3=a1+0.618*(a4-a1);
fun(x_k3,x_k,s_k,a3,n);
f3=f(x_k3);
do
{
if(f2<=f3)
{
a4=a3;
a3=a2;
f3=f2;
a2=a1+0.382*(a4-a1);
fun(x_k2,x_k,s_k,a2,n);
f2=f(x_k2);
}
else
{
a1=a2;
a2=a3;
f2=f3;
a3=a1+0.618*(a4-a1);
fun(x_k3,x_k,s_k,a3,n);
共轭梯度法C语言实现
/*------------更新 r ----------------*/
for(r2=0.0,i=0;i<n;i++){
for(rA[i]=0.0,j=0;j<n;j++)
rA[i]+=A[i][j]*x[j]; r[i]=b[i]-rA[i];
for(j=0;j<n;j++) fin>>A[i][j]; } for(j=0;j<n;j++) fin>>b[j]; fin.close();
/*--------------------------------*/
第1页
共轭梯度法.txt printf("迭代次数 = %d\n",Gong(A,b,x,n)); print(x,n);
for(int i=0;i<n;i++){ for(int j=0;j<n;j++)
printf("%f\t",A[i][j]); printf("\n");
} }
第5页
/*--------------------------------------*/
}
/*-----------计算步长 a ----------------*/
for(r2=0.0,r1=0.0,i=0;i<n;i++){
for(rA[i]=0.0,j=0;j<n;j++) 第3页
共轭梯度法.txt
for(i=0;i<n;i++)x[i]=0.0;
共轭梯度法求极小值
应用共轭梯度法求解方程组121242x x x x +=⎧⎨-=⎩的根。
初始值(0)(0,0)T x = 分析：将方程组Ax b =转化为优化问题中的极值问题然后应用共轭梯度法进行求解。
令()R x Ax b =-，只需求得x ，使得()R x 取得最小值。
则有： min ()R x ⇔2min ()R x min T R R ⇔()()()()()()()()2T T T T T T T T T T T T T T T T T T T T T T T T T T R R Ax b Ax b X A b Ax b X A AX X A b b Ax b bX A AX X A b b Ax b bX A A X b AX b Ax b bX A A X b AX b b=--=--=--+=--+=--+=-+这是一个数 与1()2T T f x x Gx b x c =++比较 则有：2,()22T T T G A A g x A A A b ==- 回到题目问题中，则对应有124124012,,04444x G g b x --⎛⎫⎛⎫⎛⎫=== ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭也就是应用共轭梯度法求点使221212()2212420f x x x x x =+--+取得最小值。
程序清单：#include <stdio.h>#include <math.h>double a,b,s[2]; /*全局变量*/double E=1e-6;double F(double x1,double x2){double y;y=2*x1*x1+2*x2*x2-12*x1-4*x2+20;return (y);}void qujian(double x1,double x2) /*用前进后退法求a 的探索区间*/ {double a0=0,h=1,a1,a2,f1,f2,X1,X2,X3,X4;a1=a0;a2=a0+h;X1=x1+a1*s[0];X2=x2+a1*s[1];f1=F(X1,X2);X3=x1+a2*s[0];X4=x2+a2*s[1];f2=F(X3,X4);if (f1>f2){while(1){h=h*2;a2=a2+h;f1=f2;X3=x1+a2*s[0];X4=x2+a2*s[1];f2=F(X3,X4);if (f1>f2) {a1=a2-h;}else{a=a1;b=a2;break;}}}else{h=-h/4;while(1){a1=a1+h;f2=f1;X1=x1+a1*s[0];X2=x2+a1*s[1];f1=F(X1,X2);if(f1<f2) {a2=a1-h;h=2*h;}else{a=a1;b=a2;break;}}}}double MIN(double x1,double x2) /*二次插值法求a的最小值*/ {double x01,x02,x03,x0,xmin;double f1,f2,f3,f;double X1,X2,X3,X4,X5,X6,X7,X8;double h=(b-a)/2;double k1,k2;x01=a;x02=a+h;x03=b;X1=x1+x01*s[0];X2=x2+x01*s[1];X3=x1+x02*s[0];X4=x2+x02*s[1];X5=x1+x03*s[0];X6=x2+x03*s[1];f1=F(X1,X2);f2=F(X3,X4);f3=F(X5,X6);while (1){k1=(f3-f1)/(x03-x01);k2=((f2-f1)/(x02-x01)-k1)/(x02-x03);x0=(x01+x03-k1/k2)/2;X7=x1+x0*s[0];X8=x2+x0*s[1];f=F(X7,X8);if(fabs(x0-x02)<=E){xmin=x0;break;}else if(x02<x0){if(f2<f){x03=x0;f3=f;}else{x01=x02;f1=f2;x02=x0;f2=f;}}else{if(f2<f){x01=x0;f1=f;}else{x03=x02;f3=f2;x02=x0;f2=f;}}}return(xmin);}void main(){double g[3],m1,m2,b;int n=2,k=0;double x1=0,x2=0; /*起始点*/g[1]=4*x1-12;g[2]=4*x2-4;m1=g[1]*g[1]+g[2]*g[2];s[0]=-g[1];s[1]=-g[2];while(1){qujian(x1,x2);a=MIN(x1,x2); /*求a的探索区间*/x1=x1+a*s[0];x2=x2+a*s[1];g[1]=4*x1-12;m2=g[1]*g[1]+g[2]*g[2];if(sqrt(g[1]*g[1]+g[2]*g[2])<=E)break;if(k+1==n){g[1]=4*x1-12;g[2]=4*x2-4;s[0]=-g[1];s[1]=-g[2];}else{b=m2/m1;s[0]=-g[1]+b*s[0];s[1]=-g[2]+b*s[1];k=k+1;}printf("最优解：\n X1=%f\n X2=%f\n",x1,x2);}}运行结果：最优解：X1=3.000000X2=1.000000。
共轭梯度法求三元函数极小值
/*//----------------以下是测试部分
Vector v1(1,2,0),v2(1,2,3),v3(3,1,2),v4,v5,v6,v7(0,3,4);double n;
Vector v_arry[133]={v1,v2,v3};
}
double GetModulus(Vector &p)
{ቤተ መጻሕፍቲ ባይዱ
return sqrt(p.x1*p.x1+p.x2*p.x2+p.x3*p.x3);
}
/*double GetFValue(double la) //求函数值f(x[k]+lamda*s[k])
{
return ValueFunction(x[k]+la*s[k]);
la0=la1;
h=h*2;
goto step12;
}
else
{
//cout<<"f1>f0"<<endl;
if(fabs(h)<e)
{
lamdak=la0;
//cout<<"|h|为："<<fabs(h)<<endl;
//cout<<"e为："<<e<<endl;
//cout<<"达到结束条件：|h|<e，计算结束。"<<endl;
double f1,f0,la1,la0=1,h=10,e=0.0000001;
共轭梯度法c++程序
共轭梯度法c++程序共轭梯度法程序源代码 #include<stdio.h> #include<math.h> #define N 10 #define eps pow(10,-6)double f(double x[],double p[],double t){double s;s=pow(x[0]+t*p[0],2)+25*pow(x[1]+t*p[1],2);return s;}/*以下是进退法搜索区间源程序*/ void sb(double *a,double *b,doublex[],double p[]){double t0,t1,t,h,alpha,f0,f1;int k=0;t0=2.5; /*初始值*/h=1; /*初始步长*/alpha=2; /*加步系数*/f0=f(x,p,t0);t1=t0+h;f1=f(x,p,t1);while(1){if(f1<f0){h=alpha*h; t=t0;t0=t1; f0=f1;k++;}else{if(k==0){h=-h;t=t1;}else{*a=t<t1?t:t1;*b=t>t1?t:t1;break;}}t1=t0+h;f1=f(x,p,t1);}}/*以下是黄金分割法程序源代码*/ double hjfg(double x[],double p[]) {double beta,t1,t2,t; double f1,f2;double a=0,b=0;double *c,*d;c=&a,d=&b;sb(c,d,x,p);/*调用进退法搜索区间*/printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf",x[0],x[1],p[0],p[1]); printf("\n[a,b]=[%lf,%lf]",a,b);beta=(sqrt(5)-1.0)/2; t2=a+beta*(b-a); f2=f(x,p,t2);t1=a+b-t2; f1=f(x,p,t1);while(1){if(fabs(t1-t2)<eps)break;else{if(f1<f2){t=(t1+t2)/2;b=t2; t2=t1;f2=f1; t1=a+b-t2;f1=f(x,p,t1);}else{a=t1; t1=t2;f1=f2;t2=a+beta*(b-a);f2=f(x,p,t2);}}}t=(t1+t2)/2;return t;}/*以下是共轭梯度法程序源代码*/void gtd(){double x[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0; int i,k,n;printf("请输入函数的元数值n=");scanf("%d",&n);printf("\n请输入初始值:\n");for(i=0;i<n;i++)scanf("%lf",&x[i]);f0=f(x,g,t);g[0]=2*x[0]; g[1]=50*x[1];mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/ if(mod1>eps){p[0]=-g[0]; p[1]=-g[1]; k=0;while(1){t=hjfg(x,p);/*调用黄金分割法求t的值*/printf("\np1=%lf,p2=%lf,t=%lf",p[0],p[1],t);x[0]=x[0]+t*p[0]; x[1]=x[1]+t*p[1];g[0]=2*x[0]; g[1]=50*x[1];/*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf",x[0],x[1],g[0],g[1]);*/ mod2=sqrt(pow(g[0],2)+pow(g[1],2)); /*求梯度的长度*/if(mod2<=eps) break;else{if(k+1==n){g[0]=2*x[0]; g[1]=50*x[1];p[0]=-g[0]; p[1]=-g[1]; k=0;}else{nanda=pow(mod2,2)/pow(mod1,2);printf("\nnanda=%lf,mod=%lf",nanda,mod2);p[0]=-g[0]+nanda*p[0];p[1]=-g[1]+nanda*p[1];mod1=mod2;k++;}}printf("\n--------------------------");}}printf("\n最优解为x1=%lf,x2=%lf",x[0],x[1]); printf("\n最终的函数值为%lf",f(x,g,t));}main(){gtd();}运行结果如下:请输入函数的元数值n=2请输入初始值:2 2x1=2.000000,x2=2.000000,p1=-4.000000,p2=-100.000000[a,b]=[-4.500000,1.500000]p1=-4.000000,p2=-100.000000,t=0.020030 nanda=0.001474,mod=3.842730 --------------------------x1=1.919879,x2=-0.003022,p1=-3.845655,p2=0.003665[a,b]=[-4.500000,1.500000]p1=-3.845655,p2=0.003665,t=0.499240 --------------------------x1=-0.000026,x2=-0.001192,p1=0.000052,p2=0.059610[a,b]=[-4.500000,1.500000]p1=0.000052,p2=0.059610,t=0.020000 nanda=0.000000,mod=0.000050 --------------------------x1=-0.000025,x2=-0.000000,p1=0.000050,p2=0.000001[a,b]=[-4.500000,1.500000]p1=0.000050,p2=0.000001,t=0.495505 -------------------------- x1=-0.000000,x2=0.000000,p1=0.000000,p2=-0.000023[a,b]=[-4.500000,1.500000]p1=0.000000,p2=-0.000023,t=0.020007 最优解为x1=-0.000000,x2=-0.000000 最终的函数值为0.000000Press any key to continue。
FR共轭梯度法实验及其程序代码
实验报告2FR共轭梯度法1 实验目的掌握外法函数法2 实验内容外罚函数法3 算法设计①编写主函数：function [min,answ]=OuterPenaltyFunction(x0,eps)实现外罚函数法②定义function y=PFun(x)用以实现求目标函数的一阶导数；③定义测试函数: function y=Fun(x)。
4 程序代码（一）外罚函数法程序代码function [min,answ]=OuterPenaltyFunction(x0,eps)%% ---- 该程序用以实现外罚函数法%% Input：% x0 ---- 初始点% eps ---- 搜索精度%% output：% min ---- 目标函数最优值% answ ---- 最优搜索点syms x1 x2xk=x0;c=10; %罚因子放大系数33deta=0.1; %罚因子while (1)fun=Fun([x1 x2])+deta*PFun([x1,x2]); %惩罚函数f=inline(fun);[min,answ]=FR_Conjugate_gradient(f, xk,eps);%调用FR迭代法函数求确定罚函数最优解if double(deta*PFun(answ))<epsbreakenddeta=c*deta;xk=answ;endansw=xk;min=double(Fun(xk));function y=PFun(x)%% ---- 该函数用以定义惩罚项函数%% Input：% x ---- 自变量%% output：% y ---- 因变量syms x1 x2syms x1 x2x1=x(1);% x2=x(2);y=(x1+1)^2;function y=Fun(x)%% ---- 该函数用以定义测试函数%% Input：% x ---- 自变量%% output：% y ---- 因变量syms x1 x2x1=x(1);x2=x(2);y=x1^2+x2^2;（二）FR共轭梯度法程序代码：function [min,answ]=FR_Conjugate_gradient(Fun, x0,eps)%% ---- 改程序用以实现FR共轭梯度法%% Input：% x0 ---- 初始点% eps ---- 搜索精度%% output：% min ---- 目标函数最优值% answ ---- 最优搜索点syms beta alphk=1;%迭代次数xk=x0;temp=1;while(norm(DFun(Fun,xk))>eps)if (k==1)pk=-DFun(Fun,xk);elsebeta=DFun(Fun,xk)'*DFun(Fun,xk)/temp;pk=-DFun(Fun,xk)+beta*pk;endtemp=DFun(Fun,xk)'*DFun(Fun,xk);temp2=xk+alph*pk';fun=Fun(temp2(1),temp2(2));fun=inline(fun);[min,answ]=OneDimentionSearch(fun,0,0.08,eps);%调用精确一维搜索求步长xk=xk+answ*pk';%answ即为所求的步长k=k+1 ;endansw=xk;min=double(Fun(xk(1),xk(2)));function y=DFun(Fun,x)%% ---- 改函数用以实现求目标函数的一阶导数%% Input：% x ---- 求导点%% output：% y ---- 目标函数的一阶导数在x点的值syms x1 x2f=Fun(x1,x2);f1=diff(f,x1);f2=diff(f,x2);f1=subs(f1,{x1,x2},x);f2=subs(f2,{x1,x2},x);y=[double(f1),double(f2)]';（三）精确一维搜索程序代码：function [min,answ]=OneDimentionSearch(Fun,x0,step,eps) %% ---- 改程序用以实现一维搜索算法%% Input：% Fun ---- 一维精确搜索函数% x0 ---- 初始点% step ---- 搜索步长% eps ---- 搜索精度%% output：% min ---- 目标函数最优值% answ ---- 最优搜索点%%% 求初始区间[a,b]=AdvanceAndRetreat(Fun,x0,step);%黄金分割法实现一维搜索[min,answ] =GoldSection(Fun,a,b,eps);function [min,answ] =GoldSection(Fun,a,b,eps)%% ---- 改函数用以实现黄金搜索法%% Input：% Fun ---- 一维精确搜索函数% a ---- 搜索区间的左端点% b ---- 搜索间隔的右端点% eps ---- 搜索精度%% output：% min ---- 目标函数最优值% wnsw ---- 最优搜索点%%format longx1=a+.382*(b-a);x2=a+.618*(b-a);while(abs(b-a)>eps)f1=Fun(x1);f2=Fun(x2);% 比较判断两个分割点处的函数值，进而缩短区间长度if(f1>f2)a=x1; x1=x2;x2=a+.618*(b-a);elseif(f1==f2)a=x1; b=x2;x1=a+.382*(b-a);x2=a+.618*(b-a);elseb=x2;x2=x1;x1=a+.382*(b-a);endend% 返回搜索点和搜索值answ=(a+b)/2;min=Fun(answ);function [a,b]=AdvanceAndRetreat(Fun,x0,step)%% ---- 改函数用以实现进退法求初始区间%% Input：% Fun ---- 一维精确搜索函数% x0 ---- 初始点% step ---- 搜索间隔%% output：% a ---- 搜索区间的左端点% b ---- 搜索间隔的右端点x1=x0+step;if Fun(x1)<=Fun(x0)while(1)step=2*step;x2=x1+step;if Fun(x1)<=Fun(x2)break;else x0=x1;x1=x2;endenda=x0;b=x2;elsewhile(1)step=2*step;x2=x0-step;if Fun(x0)<=Fun(x2)break;else x1=x0;x0=x2;endenda=x2;b=x1;end5 运行结果测试函数为：01..)(min 12221=++=x t s x x x f该函数是教材第147页例4.2.2，其结果与书上的结果近似相等（因为存在微小迭代误差）。