• ## JACOBI

2019-10-24 12:51:03
JACOBI
• ## Jacobi

2018-11-05 13:55:29
Jacobi function[x,k,index]=Jacobi(A,b,ep,it_max) if nargin&lt;4 it_max=100; end if nargin&lt;3 ep=1e-5; end n=length(A); k=0; x=zeros(n,1); y=zeros(n,1); index=1; while k&lt;=it_max ...
Jacobi

function[x,k,index]=Jacobi(A,b,ep,it_max)
if nargin<4
it_max=100;
end
if nargin<3
ep=1e-5;
end

n=length(A);
k=0;
x=zeros(n,1);
y=zeros(n,1);
index=1;
while k<=it_max
for i =1:n
if abs(A(i,i))<1e-10
index=0;
return;
end
y(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i))/A(i,i);
end
if (x,inf)<ep
break;
end
k=k+1;
x=y;
end


展开全文
• ## jacobi

2015-05-11 23:28:08
SUBROUTINE jacobi(a,n,np,d,v,nrot) ! ! Computes all eigenvalues and eigenvectors of a real symmetric matrix a, ! which is of size n by n, stored in a physical np by np array. ! On output, elements o
      SUBROUTINE jacobi(a,n,np,d,v,nrot)
!
! Computes all eigenvalues and eigenvectors of a real symmetric matrix a,
! which is of size n by n, stored in a physical np by np array.
! On output, elements of a above the diagonal are destroyed.
! d returns the eigenvalues of a in its first n elements.
! v is a matrix with the same logical and physical dimensions as a,
! whose columns contain, on output, the normalized eigenvectors of a.
! nrot returns the number of Jacobi rotations that were required.
! Please notice that the eigenvalues are not ordered on output.
! If the sorting is desired, the addintioal routine "eigsrt"
! can be invoked to reorder the output of jacobi.
!
INTEGER n,np,nrot,NMAX
REAL a(np,np),d(np),v(np,np)
PARAMETER (NMAX=500)
INTEGER i,ip,iq,j
REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do 12 ip=1,n
do 11 iq=1,n
v(ip,iq)=0.
11      continue
v(ip,ip)=1.
12    continue
do 13 ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.
13    continue
nrot=0
do 24 i=1,50
sm=0.
do 15 ip=1,n-1
do 14 iq=ip+1,n
sm=sm+abs(a(ip,iq))
14        continue
15      continue
if(sm.eq.0.)return
if(i.lt.4)then
tresh=0.2*sm/n**2
else
tresh=0.
endif
do 22 ip=1,n-1
do 21 iq=ip+1,n
g=100.*abs(a(ip,iq))
if((i.gt.4).and.(abs(d(ip))+
*g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
a(ip,iq)=0.
else if(abs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(abs(h)+g.eq.abs(h))then
t=a(ip,iq)/h
else
theta=0.5*h/a(ip,iq)
t=1./(abs(theta)+sqrt(1.+theta**2))
if(theta.lt.0.)t=-t
endif
c=1./sqrt(1+t**2)
s=t*c
tau=s/(1.+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.
do 16 j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
16            continue
do 17 j=ip+1,iq-1
g=a(ip,j)
h=a(j,iq)
a(ip,j)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
17            continue
do 18 j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(ip,j)=g-s*(h+g*tau)
a(iq,j)=h+s*(g-h*tau)
18            continue
do 19 j=1,n
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
19            continue
nrot=nrot+1
endif
21        continue
22      continue
do 23 ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.
23      continue
24    continue
pause 'too many iterations in jacobi'
return
END
</pre><pre name="code" class="cpp">

展开全文
• Matlab计算Jacobi、Gauss-Seidel迭代例题Jacobi迭代Gauss-Seidel迭代 例题 方程组{5x1+2x2+x3=−12−x1+4x2+2x3=202x1−3x2+10x3=3\begin{cases} 5x_1 + 2x_2 + x_3 = -12\\ -x_1 + 4x_2 + 2x_3 = 20\\ 2x_1 - 3x_2 ...
使用Matlab实现：Jacobi、Gauss-Seidel迭代例题Jacobi迭代Gauss-Seidel迭代
例题
方程组$\begin{cases} 5x_1 + 2x_2 + x_3 = -12\\ -x_1 + 4x_2 + 2x_3 = 20\\ 2x_1 - 3x_2 + 10x_3 = 3\\ \end{cases}$	求解，当$max|x_i^{(k + 1)} - x_i^{(k)}| \leq 10^{-5}$ 时候迭代终止。
以下解答过程，上标表示迭代次数，下标表示序号。
Jacobi迭代
定义变量：
$D = diag(a_{11}, a_{22}, ..., a_{nn}),\\ L = \left[\begin{array}{cccccc} 0\\ -a_{21} & 0\\ ...\\ -a_{i1} & ... & -a_{i,i-1} & 0\\ ...\\ -a_{n1} & ... & -a_{n,i-1} & ... & -a_{n,n-1} & 0\\ \end{array}\right],\\ U = \left[\begin{array}{cccccc} 0 &-a_{12} & ... & -a_{1,i} & ... & -a_{1,n}\\ ...\\ & & 0 & -a_{i-1,i} & ... & -a_{i-1,n}\\ ...\\ & & & & 0 & -a_{n-1,n}\\ ...\\ & & & & & 0\\ \end{array}\right]$
其矩阵迭代形式为：
$x^{(k+1)} = B_J \cdot x^{(k)} + f_J\\ B_J = D^{-1} \cdot (L + U), \quad f_J = D^{-1} \cdot b$。
写出分量形式：
$\begin{cases} x_1^{(k + 1)} = \frac15 (-12 - 2 x_2^{(k)} - x_3^{(k)} )\\ x_2^{(k + 1)} = \frac14 (20 + x_1^{(k)} - 2x_3^{(k)} )\\ x_3^{(k + 1)} = \frac1{10} (3 - 2 x_1^{(k)} + 3x_2^{(k)} )\\ \end{cases}$
写成矩阵形式：
$\left[\begin{array}{c} x_1^{(k + 1)}\\ x_2^{(k + 1)}\\ x_3^{(k + 1)}\\ \end{array}\right] = \left[\begin{array}{cccc} 0 & -0.4 & -0.2\\ 0.25 & 0 & -0.5\\ -0.2 & 0.3 & 0\\ \end{array}\right] \cdot \left[\begin{array}{c} x_1^{(k)}\\ x_2^{(k)}\\ x_3^{(k)}\\ \end{array}\right] + \left[\begin{array}{c} -2.4\\ 5\\ 0.3\\ \end{array}\right]$
取初始向量：$x^0 = (0, 0, 0)^T$，依次按照上式进行迭代。使用Matlab进行编程求解。
a=[0,-0.4,-0.2;0.25,0,-0.5;-0.2,0.3,0];
b = [-2.4;5;0.3];
x = [0;0;0];
xx = a * x + b;
i = 0;
while norm(x - xx, inf) >= 1e-5
x = xx;
xx = a * x + b;
i = i +1;
end

以上代码，最终 $x = x^{i}, xx = x^{(i + 1)}$ ，最终迭代次数位 $i + 1$ 次，如果你需要看到更长的小数位置，可以使用以下Matlab代码，表示使用15位浮点或定点数。
format long g

运行结果为：

即精确解为 $x = (-4,3,2)^T$ 。
Gauss-Seidel迭代
其矩阵迭代形式为：
$x^{(k+1)} = B_G \cdot x^{(k)} + f_G\\ B_G = (D - L) ^{-1} \cdot U, \quad f_G = (D - L) ^{-1} \cdot b$
使用Matlab编程求解：
d = [5,0,0;0,4,0;0,0,10];
l = [0,0,0;1,0,0;-2,3,0];
u = [0,-2,-1;0,0,-2;0,0,0];
b = [-12;20;3];
t = inv(d - l);
bg = t * u;
fg = t * b;
x = [0;0;0];
xx = [-2.4;4.4;2.1];
i = 0;
while norm(x - xx, inf) >= 1e-5
x = xx;
xx = bg * x + fg;
i = i +1;
end

运行结果为：

同样求得精确解为 $x = (-4,3,2)^T$ 。


展开全文
• 本科课程参见：《软件学院那些课》 牛顿迭代公式 设已知方程f(x)=0的近似根x0 ,则在x0附近f(x)可用一阶泰勒多项式近似代替....用x1表示p(x)=0的根,它与f(x)=0的根差异不大. 设 ,由于x1满足解得 ...
本科课程参见：《软件学院那些课》

牛顿迭代公式

设已知方程f(x)=0的近似根x0 ,则在x0附近f(x)可用一阶泰勒多项式近似代替.因此, 方程f(x)=0可近似地表示为p(x)=0。用x1表示p(x)=0的根,它与f(x)=0的根差异不大.

设 ,由于x1满足解得

重复这一过程,得到迭代公式：

这就是著名的牛顿迭代公式,它相应的不动点方程为

Jacobi迭代公式解线性方程组

线性方程组基本解法：

若方程组可同解变形为

Jacobi迭代法的计算公式：

即

算法代码

/*简单迭代法的代码实现*/
#include<iostream>
#include<string>
#include<cmath>
using namespace std;
double e=2.718281818284;
double f(double x){
return pow(e,-1*x);
}
void SimpleDiedai(double x,double d){
double a=x;
double b=f(a);
int k=0;//记录循环的次数
while(((a-b)>d) || ((a-b)<-1*d)){
cout<<a<<endl;
a=b;
b=f(a);
k++;
if(k>100){
cout<<"迭代失败！（可能是函数不收敛）"<<endl;
return ;
}
}
cout<<b<<endl;
return;
}
int main(){
cout<<"请输入初始值x0和要求得结果的精度：";
double x,d;
cin>>x>>d;
SimpleDiedai (x,d);
return 0;
}

/*牛顿迭代法的代码实现*/
#include<iostream>
#include<string>
#include<cmath>
using namespace std;
double e=2.718281818284;
double f(double x){
double a=pow(e,-1*x);
return x-(x-a)/(1+a);
}
void NewtonDiedai(double x,double d){
double a=x;
double b=f(a);
int k=0; //记录循环的次数
while(((a-b)>d) || ((a-b)<-1*d)){
cout<<a<<endl;
a=b;
b=f(a);
k++;
if(k>100){
cout<<"迭代失败！（可能是函数不收敛）"<<endl;
return ;
}
}
cout<<b<<endl;
return;
}

int main(){
cout<<"请输入初始值x0和要求得结果的精度：";
double x,d;
cin>>x>>d;
NewtonDiedai(x,d);
return 0;
}

/*雅可比算法的代码实现*/
#include<iostream>
#include<iomanip>
#include<string>
#include<vector>
using namespace std;

//函数求数组中的最大值
double MaxOfList(vector<double>x){
double max=x[0];
int n=x.size();
for(int i=0;i<n;i++)
if(x[i]>max) max=x[i];
return max;
}

//雅可比迭代公式
void Jacobi(vector<vector<double> > A,vector<double> B,int n){
vector<double> X(n,0);
vector<double> Y(n,0);
vector<double> D(n,0);
int k=0; //记录循环次数
do{
X=Y;
for(int i=0;i<n;i++){
double tem=0;
for(int j=0;j<n;j++){
if(i!=j) tem += A[i][j]*X[j];
}
Y[i]=(B[i]-tem)/A[i][i];
cout<<left<<setw(8)<<Y[i]<<" ";
}
cout<<endl;
k++;
if(k>100){
cout<<"迭代失败！（可能是函数不收敛）"<<endl;
return ;
}

for(int a=0;a<n;a++){
D[a]=X[a]-Y[a];
}
}while( MaxOfList(D)>0.00001 || MaxOfList(D)<-0.00001);

return ;
}

int main(){

int n;
cout<<"请输入方程组未知数的个数n：";
cin>>n;
cout<<endl;

vector<vector<double> >A(n,vector<double>(n,0));
vector<double>B(n,0);

cout<<"请输入方程组的系数矩阵："<<endl;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
cin>>A[i][j];
}
}
cout<<endl;

cout<<"请输入方程组的值向量："<<endl;
for(int k=0;k<n;k++){
cin>>B[k];
}
cout<<endl;

cout<<"您输入的方程组为："<<endl;
for(int a=0;a<n;a++){
for(int b=0;b<n;b++){
cout<<A[a][b]<<" ";
}
cout<<"    "<<B[a]<<endl;
}
cout<<endl;
cout<<"由雅可比迭代公式求的方程组的解为："<<endl;
Jacobi(A,B,n);
return 0;
}


实验过程原始记录

（1）分别用简单迭代法和牛顿迭代法求方程在x=0.5附近的一个根x*，要求精度为0.00001

（输入0.5 0.000001）简单迭代法得到结果：

（输入0.5 0.000001）牛顿迭代法得到结果：
X0=0.5  x1=0.566311 x2=0.567143

（2）用雅可比迭代法求解方程组

运行程序，根据提示输入 （3） （10 -1 -2 -1 10 -2 -1 -2 5）    （7.2  8.3  4.2）

实验结果及分析

1、迭代法是一种逐次逼近法，这种方法使用某个固定公式－所谓迭代公式反复校正根的近似值，使之逐步精确化，直至满足精度要求的结果。迭代法是一种求解函数方程f（x）=0的解的方法，他解决了二分法无法求解复根级偶重根的问题，而其提高了收敛速度。迭代的思想是计算方法中基础的求解问题的思想。
2、简单迭代法的求根过程分成两步，第一步先提供根的某个猜测值，即所谓迭代初值，然后将迭代初值逐步加工成满足精度要求的根。迭代法的设计思想是：f (x) = 0等价变换成 然后由迭代公式 逐步球的满足精度的解。实际迭代中不同迭代函数的求解可能影响求的精确解的运算量，甚至可能因为函数发散而无法求解。解题时可通过对导函数的判断而判断函数是否发散，而编写代码时可以通过判断循环次数——即循环过多次而不能从循环中出来时就判断为死循环，无法求得正解
3、简单迭代法往往只是线性收敛，为得出超线性收敛的迭代格式，通常采用近似替代法， 即牛顿公式。迭代函数为  - / 牛顿法是一种逐步线性化方法。由实验结果可以看到，虽然选取近似公式，但牛顿迭代法仍能得到精度很高的解，而且牛顿迭代法大大提高了收敛速度。
4、由迭代法求解线性方程组的基本思想是将联立方程组的求解归结为重复计算一组彼此独立的线性表达式，这就使问题得到了简化，类似简单迭代法转换方程组中每个方程式可得到雅可比迭代式

迭代法求解方程组有一定的局限性，例如要求方程组的系数矩阵具有某种特殊性质，以保证迭代过程的收敛性，但迭代法同时有十分明显的优点——算法简单，因而编制程序比较容易，所以在实际求解问题中仍有非常大利用价值。

（转载请注明作者和出处：http://blog.csdn.net/xiaowei_cqu 未经允许请勿用于商业用途）


展开全文
• 数值分析1，即数值代数所用程序，用Jacobi迭代法求解方程足
• 使用Jacobi–Lie群和Jacobi–Lie双代数的概念，我们将泊松–李对称性的定义推广到Jacobi–Lie对称性。 在这方面，我们将泊松－李T对偶的概念推广到雅可比－李T对偶，并给出了具有雅可比－李对称性的李群上的雅可比－...
• 为此，我们使用Jacobi–Maupertuis Randers–Finsler度量将广义光学度量方法扩展到广义Jacobi度量方法。 更具体地说，我们将Gauss-Bonnet定理应用于广义Jacobi度量空间，然后获得一个用于计算偏转角的表达式，该...
• ## Jacobi迭代法与Gauss-Seidel迭代法

万次阅读 多人点赞 2016-01-25 19:29:58
今天刚好有朋友和我讨论泊松图像融合算法，我说我过去文章里给出的是最原始、最直观的实现算法。对于理解泊松融合的原理比较有帮助，但是效率可能并不理想。印象中，泊松融合是有一个以矩阵为基础的快速算法的。...
• matlab代码，Jacobi迭代算法，亲测可用，直接调用该函数即可。
• <p>Replace the jacobi() method in math_extra with a new verison(s). The old one was an altered version of a Numerical Recipes method. This one is not and can thus be released GPL or LGPL. <p><strong>...
• hamilton jacobi方程的详细数学解释和介绍。
• 该方法主要利用移位的Jacobi多项式将方程中的函数逼近,再结合Captuo类型的变分数阶微积分定义,推导出移位Jacobi多项式的微积分算子矩阵,将最初的方程转化为矩阵相乘的形式,然后通过离散变量,将原方程转化为一系列非...
• 进过18次Jacobi迭代后，其相邻迭代解间无穷范数误差小于：1.0e-8 此时Jacobi迭代解如下： x = 1.099999996412137 1.199999996412137 1.299999995744652
• <ul><li>Less confusing variable naming</li><li>Use random Z coordinates for the (initial) gej variables (instead of 1).</li><li>Vary the inputs to the Jacobi symbol benchmarks (they were constantly ...
• <p>I will update this PR with an implementation of the Jacobi symbol computation; for now I'm pushing it in the hopes that somebody will see some more perf improvements that will let us catch up ...
• % Jacobi方法和Gauss-Seidel方法 clear all A = [8,-1,1;2,10,-1; 1 1 -5]; b = [1;4;3]; n = length(b); itermax = 100; iter = 0; tol = 1e-4 x = zeros(n,1); while iter  iter = iter + 1;  for ...
• Hamilton Jacobi 使用的不同的 Flux 和不同的边界条件测试了 $u_t=\frac{u_x^2}{1+u_x^2}-\frac{\cos(x)^2}{1+\cos(x)^2},x\in [0,4\pi]$ 最有趣的是测试了初值 \[ u= \left\{ \begin{array}{c} \sin(x),&...
• 求微分其实就是线性化，导数其实就是线性空间之间的线性变换，Jaocibian矩阵本质上就是导数。 比如，映射在处的导数就是在处的切空间到在处的切空间之间的线性映射。切空间都是矢量空间，都有基底，所以这个线性...
• ## Jacobi 方法

千次阅读 2013-11-05 10:59:03
第三节 Jacobi 方法  Jacobi方法是求对称矩阵的全部特征值以及相应的特征向量的一种方法，它是基于以下两个结论  1) 任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q,使得  QT AQ = diag(λ1 ...
• #include<stdio.h> #include<math.h> #define N 3 void print(double s[N][N]) { for (int i = 0; i <N; i++) { for (int j = 0;...void multiply(double g[N][N],double v[N][N])//乘
• (1)Jacobi 迭代法。 (2)Gauss-Seidel 迭代法。 (3)逐次超松弛迭代法。 (4)共轭梯度法。 A 为对称正定矩阵，其特征值服从独⽴同分布的[0,1]间的均匀分布;b 中的元素服从独立同 分布的正态分布。令 n=10、50、100、200...
• Jacobi 迭代法C语言实现，已上机调试无误，三阶方程组，若需求解高阶方程组，可根据提示小改即可。
• <p>As this leaves the entire num/gmp support only required for the (unused) Jacobi symbols, remove support for all of those entirely.</p><p>该提问来源于开源项目：bitcoin-core/secp256k1</p></div>
• 并行计算课程作业，实现jacobi迭代的串行优化，主要是一级和二级cache优化。
• ## Jacobi迭代法

万次阅读 2018-08-08 17:02:15
其中Jacobi迭代法就是很简单易懂的一种。 我编写的C++代码如下，其中文件matrix.txt内容如下， 第一行的第一个数字表示矩阵的行数或者列数，因为rows==cols 下面的三行表示矩阵本体 最后的一行表示该矩阵和向量...

...