Can someone help me setup the Chandrasekhar H-equation in Matlab?

I have a couple errors I cannot seem to resolve. Below is the discrete equation form.

This is a homework problem so I am not looking for answers, just help to setup the problem which I will solve using various methods.

First I setup the internal structure

global ChH;
c=.8;
N=100;

meu=1:N;
meu=(meu-.5)/N;
meu=meu';

cs=.5*c/N;
ChH=ones(N,1)*meu'; 

ChH=cs*ChH'./(ChH+ChH');

Then I setup the rest of the function:

   function [ x, ex ] = Heq_Nwt( f, df, x1, tol, nmax )   

global ChH;

N=100;
n2=ones(N,N);
x1 = ones(N,N);

ch2=n2-(ChH.*(x1'));      % [1 - c/2N (sum: meu*x_j /(meu+meu'))]
sm= n2./ch2;            % 1/previous
ch2=x1-sm                % x - previous

Problem is I am not sure what my starting value of x should be? ChH is 100x100 so x must be same dimensions for multiplication right? I converted 1 to a 100x100 matrix as well for subtraction and division? I did not include any sort of summation for internal N_summation_ j=1? Shouldn't x be a vector of length N per the equation statement given? Then why is equation structure creating matrix form?

Answer 1

[Ahmet Cecen]

I will tell you 2 things:

1) x is most definitely a vector.

2) The NxN matrix is obtained by writing the summation over j as a matrix vector product, and a matrix vector product is itself a vector.

[Branden Chamness]

How can I setup x unknown vector to use in equation through summation? Use symsum as below?

sym x
c=.8;
N=100;
meu=1:N;
meu=(meu-.5)/N;
meu=meu';
cs=.5*c/N;
ChH= meu';

symsum(x, 0, N)

ChH=symsum(x*(ChH'./(ChH+ChH'), 0, 10));

ChH = ChH*cs

[Ahmet Cecen]

http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html

[Ahmet Cecen]

You are going about the implementation in a very weird way, I am not even sure you can solve it that way. The summation should never appear explicitly in your script, replace it with a matrix vector product. Paste following in a new script and publish:

 %%
 %
 % $$ (Ax)_i = {{c}\over{2N}}\;\sum_{j=1}^{N} {{\mu_i}\over{\mu_i+\mu_j}}x_i $$
 %
 % $$ A_{ij} = ({{c}\over{2N}}) ({{\mu_i}\over{\mu_i+\mu_j}}) $$ 

Answer 2

[Branden Chamness]

I used the following for summation and I am not receiving any errors

in = @(x1) sum(x1.*(ChH'./(ChH+ChH')))

I'd like to do something like the following to finish the function but not sure how to use above function handle to input for final equation.

f = @(x1)x1' - inv(n2-(cs*in));