1

I'm trying to use fsolve in matlab to solve a system of nonlinear equations numerically. Here is a test sample of my program, k1 and R are parameters and x0 is the start point.

function y=f(k1, R, x0)

pair=fsolve(@system,x0);

y=pair(1);

    function r=system(v)

        int1=@(x) exp(k1*x);
        int2=@(x) exp(k1*x^2)/(x^4);


        r(1)=exp(v(1))*quadl(int1,0,v(1));     
        r(2)=exp(k1*v(2))*quadl(int2,v(1),20)*k1*R;

    end

end

The strange thing is when I run this program, matlab keeps telling me that I should use .^ instead of ^ in int2=@(x) exp(k1*x^2)/(x^4). I am confused because the x in that function handle is supposed to be a scalar when it is used by quadl. Why should I have to use .^ in this case?

Also I see that a lot of the examples provided in online documentations also use .^ even though they are clearly taking power of a scalar, as in here. Can anybody help explain why?

Thanks in advance.

2 Answers 2

2

in the function int2 you have used matrix power (^) where you should use element-wise power (.^). Also, you have used matrix right division (/) where you should use element-wise division (./). This is needed, since quadl (and friends) will evaluate the integrand int2 for a whole array of x's at a time for reasons of efficiency.

So, use this:

function y = f(k1, R, x0)

    pair = fsolve(@system,x0);

    y = pair(1);

    function r = system(v)

        int1 = @(x) exp(k1*x);
        int2 = @(x) exp(k1*x.^2)./(x.^4);

        r(1) = exp(   v(1)) * quadl(int1,0,v(1));     
        r(2) = exp(k1*v(2)) * k1*R*quadl(int2,v(1),20);

    end

end

Also, have a look at quadgk or integral (if you're on newer Matlab).

By the way, I assume your real functions int1 and int2 are different functions? Because these functions are of course trivial to solve analytically...

10
  • Thank you very much Rody! Yeah the real ones are different (in fact they do not even have a closerd form). But I'm getting errors in them and I'm new to matlab, so I try to decompose it into test samples and see where the bug is :)
    – Vokram
    Oct 24, 2012 at 13:37
  • @Vokram: Elliptic integrals and Bessel functions etc. are not closed-form either, but can be evaluated much faster than a quadrature method. Have you tried them in the online Mathematica integrator? Oct 24, 2012 at 13:40
  • Thanks Rody. Apologies for my limited knowledge in mathematics... The function I am trying to evaluate is an integral of an integral (from an economic model I'm building), so my current plan is to use multiple function handles in fslove. Not sure if this fits in with your method?
    – Vokram
    Oct 24, 2012 at 14:04
  • 1
    @Vokram: Well now, that makes analytic solutions possible! For the first integral, click here. For the second integral (and third and fourth), click here. The first one is by far the most complicated. You'll have to look into the hypergeom function and how to get 2F1(a,b,c,x) out of it. Oct 25, 2012 at 3:54
  • 1
    @Vokram: Note that the first integral is also not definite (e.g., you still have to evaluate it at both S and 0 and subtract). My suspicion is that, although complicated, it will be orders of magnitude faster. Especially the second equation of course becomes easy to evaluate. Oct 25, 2012 at 3:56
2

Internally MATLAB will evaluate the function fun for necessary values of x, which is more than one and x is a vector. a and b are only used to describe the limits of integration. From the documentation

fun is a function handle. It accepts a vector x and returns a vector y, the function fun evaluated at each element of x. Limits a and b must be finite.

Hence, you must use .^ to operate on individual elements of vector x.

3
  • Thank you angainor, but can you elaborate a bit more on what the "necessary values of x" are? Are they the points where quadl attempts to calculate the function numerically and see how close the answer is to zero?
    – Vokram
    Oct 24, 2012 at 13:32
  • @Vokram Those are defined for the adaptive Lobatto quadrature used in quadl (see documentation). The method is based on evaluating the function in specific points in order to integrate it numerically.
    – angainor
    Oct 24, 2012 at 13:58
  • @Vokram The x values will be from [a, b] range, but there will be many of them.
    – angainor
    Oct 24, 2012 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.