I want to solve for $n_{0}$ for a fixed value of $n$, lets say $n=1$ $$ n= n_{0}+ \dfrac{1}{2}\int_{-1/2}^{1/2}dq\left(\dfrac{e_{q}+Un_{0}}{\hbar\omega}-1\right) $$ where $e_{q}=2[1-cos(2\pi q)] $ and $\hbar\omega=(e_{q}^{2}+2Un_{0}e_{q})^{1/2} $

Then I want to plot between $U$ and $n_{0}$. The following is the code in MATLAB.

 eq = @(q) 2*(1-cos(2*pi*q));
 hq = @(q,U,n0) ((eq(q))^2+2*U*n0*(eq(q))).^(1/2);
 y = @(q,U,n0) (((eq(q))+(U*n0))/hq(q,U,n0))-1;
 a = -0.5;
 b = 0.5;
 v = @(U,n0) quad(@(q) y(q,U,n0),a,b);
 cv =@(U,n0) n0+(0.5*v)-1;
 n0 = 0.1;
 options = optimset('Display', 'iter'); 
 n = @(U) fsolve(vv,n0,0.1);
 U = 0:0.1:20;

edit 2 code corrected edit 1.

Edit 3

But I still facing the following errors:

Error using ^ Inputs must be a scalar and a square matrix. To compute elementwise POWER, use POWER (.^) instead.

Error in @(q,n0)((eq(q))^2+2*Un0(eq(q))).^0.5

Error in @(q,n0)(((eq(q))+(U*n0))/hq(q,n0))-1

Error in @(q)y(q,n0)

Error in quad (line 72) y = f(x, varargin{:});

Error in @(n0)quad(@(q)y(q,n0),a,b)

Error in @(n0)n0+(0.5*v(n0))-1

Error in simulv1 (line 12) n = fsolve(cv(n0),n0,0.1,options);

  • 3
    $\begingroup$ welcome to StackExchange. It would help if you listed the errors and provided a description of what you have already tried and what it is that you can't seem to work out. $\endgroup$ Feb 20 '15 at 21:07
  • 1
    $\begingroup$ Welcome to SciComp! As currently posed, your question is unclear. As WolfgangBangerth points out, it's not clear what you've tried. It looks like you're relying on quad to calculate the integral, and then passing that integral as a function to fsolve. The algorithms used by fsolve require derivatives of the nonlinear equation being solved. Since quad uses adaptive quadrature algorithms, it would not surprise me if the finite difference derivative approximations used by fsolve were ill-behaved. You might have better luck using a non-adaptive quadrature instead. $\endgroup$ Feb 21 '15 at 5:19
  • $\begingroup$ Your first step would be fixing the equations in your code as per the answer to your question here. $\endgroup$
    – horchler
    Feb 23 '15 at 17:45
  • $\begingroup$ The first power in your second line needs to be eq(q).^2, not eq(q)^2. (Before you edit again: the following lines have similar problems; anytime a vector such as eq(q) appears, you need .*, ./ etc.) $\endgroup$ Feb 25 '15 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.