sormatlab


for i=2:n
    p(i,1)==0.0;
    p(1,i)=0.0;
    p(1,n+1)=0.0;
    P(n+1,1)=0.0;
    for j=2:n
        p(i,j)=1/c;
        d(i,j)= lambda1+lambda2+mhu1+mhu2;
    end
end
for i=2:n
    d(2,i)= d(2,i)-mhu2;
    d(i,2) = d(i,2)-mhu1;
    d(n,i)=d(n,i)-lambda2;
    d(i,n)=d(i,n)-lambda1;
end
relerr=1.0;
k=0;
while (relerr >1.e-7 && k<100)
		k=k+1;
		diffs=0.0;
		sump=0.0;
        for i=2:n
            for j=2:n
			y=lambda2*p(i-1,j)+mhu2*p(i+1,j) +lambda1*p(i,j-1)+mhu1*p(i,j+1);
			x=y-d(i,j)*p(i,j);
			diff=omega*x/d(i,j);
			p(i,j)=p(i,j)+diff;  
			sump=sump+p(i,j);
			diffs=diffs+diff*diff;
            end
        end
        disp([k relerr]);
		relerr=diffs/sump;
    
end

sorc

#include <iostream>

using std::cout;
using std::endl;
using std::cin;   
int main()
{
	double p[100][100],d[100][100];
	int n=7		;
	double omega, lambda1,mhu1, lambda2, mhu2,c,x,diff,diffs, sump, y,relerr;
	int i,j,k;	
	cout << "n="<<n << endl;
	c=1.0/((n-1)*(n-1));
	for (i=0; i<=n;i++)
	{ p[0][i]=0.0; p[i][0]=0.0; p[i][n]=0.0 ;
	p[n][i]=0.0;
	}
	while (true)
	{
		cout <<"\n enter omega, 0 to stop"<<endl;;
		cin >> omega;
		if (omega==0.0) return 0;
		cout <<" enter lambda1, lambda2, mhu1,mhu2"<<endl;
		cin >> lambda1 >>lambda2>>mhu1>>mhu2;
		if (omega>2.0)
			break;
	//initialize solution
	for (i=1;i<n;i++)
		for (j=1;j<n;j++)
		{p[i][j]=c ;  //initial quess
		//positive diagonal so negative off diagonal
			d[i][j]=lambda1+mhu1+lambda2+mhu2; 
		}
			
	cout <<"made p"<< endl;
	// fix up diagonal around the border
	for (i=1;i<n;i++)
	{
		d[1][i]-=mhu2;
		d[n-1][i]-=lambda2;
		d[i][1]-=mhu1;
		d[i][n-1]-=lambda1;
	}
	cout <<"made d omega " <<omega <<endl;2
	relerr=1.0;
	k=0;
	while (relerr >1.e-7)
	{
		k++;
		diffs=0.0;
			sump=0.0;
	for (i=1;i<n;i++)
		for (j=1;j<n;j++)

		{	//put in y the nearest neighbors
			y=lambda2*p[i-1][j]+mhu2*p[i+1][j] +
			   lambda1*p[i][j-1]+mhu1*p[i][j+1];
			   x=y-d[i][j]*p[i][j];
			   
			diff=omega*x/d[i][j];
		
			p[i][j]=p[i][j]+diff;
			

			sump+=p[i][j];
			diffs+=diff*diff;
					}
		relerr=diffs/sump;
		cout << " k= "<< k <<" relative err "<<relerr<<endl;
	}
	cout <<"omega "<<omega << "n " <<n<< endl;
	for (i=1;i<n;i++)
		for (j=1;j<n;j++)
			p[i][j]=p[i][j]/sump;

	
    for (i=1; i<n; i++)
	{
		cout << endl;
		for (j=1; j<n; j++)
			cout << p[i][j]<<" " ;
	}

	}	
	

	return 0;
}