Using fsolve in a loop

Hey folks I have a loop that creates 4 equations and would like to use fsolve on them. I have been using solve up to now but as that solves symbolically I'm assuming fsolve will be faster..?

Each time the loop runs it produces two vectors of equations which are (for example)

eq1 = [0.25*k12 - 0.038675*k22;-9.81*sin(0.25*k11 - 0.038675*k21 + 0.7854)]

and

eq2 = [0.53868*k12 + 0.25*k22;-9.81*sin(0.53868*k11 + 0.25*k21 + 0.7854)]

then I set eq1 equal to k1 = [k11;k12] and eq2 equal to k2 = [k21;k22] and solve for k11 k12 k21 and k22. (Note that I have used vpa to cut eq1 and eq2 down to 5 decimal places, normally these have fractions in them instead.)

Now certain values in eq1 and eq2 will change every stage of the loop so I can't use a predefined function in fsolve.

So how do I set up this fsolve to find the values k11 k12 k21 k22?

I have just been doing it the way I would have expected it to work which I thought would be

x = fsolve([eq1 - k1; eq2 - k2], x0)

but that didn't work so I read up a bit on fsolve and again tried but this time I tried to define a vector F like so

F = @(k11,k12,k21,k22)[eq1 - k1; eq2 - k2];

and then do fsolve(F, x0) as I saw the @() part used in another thread although this did not work either. (I didn't think it would though since I think the @ part works if F was defined using say x(1) and x(2) then it would be @(x).

So can anyone give me any ideas as to how I could use fsolve here?

Thanks.

EDIT: Should also say the equations are found using the commands

eq1 = feval(fun, 1, out(j,:)' + 0.25*dt*k1 + (0.25 - gamma)*dt*k2);

eq2 = feval(fun, 1, out(j,:)' + (0.25 + gamma)*dt*k1 + 0.25*dt*k2);

where out, gamma and dt are already defined, fun is my input function and my loop runs from j=1:n (n defined earlier as well). k1 and k2 are defined as above.

Answer 1

[Matt Tearle]

This seems a bit confused... why are you using feval? You should be able to just evaluate the function handle fun at the given values. But then eq1 is just a value, not a function handle, so you won't be able to use it in fsolve. How about making it one 4-D problem, with a 4-D vector of unknowns k = [k11;k12;k21;k22]. Then something like this:

eq = @(k) [fun(1,out(j,:)' + 0.25*dt*k(1) + (0.25 - gamma)*dt*k(2));...
           fun(1,out(j,:)' + (0.25 + gamma)*dt*k(3) + 0.25*dt*k(4))] - k;

Apologies if that's not quite right -- I'm trying to interpret what you have. I assume your feval statements were returning 2-by-1 vectors, so I'm vertically concatenating into one 4-by-1.

Then, of course, do kout = fsolve(eq,k0)

[Brian]

I don't think I'm gonna be able to explain this one without posting my whole code and I need to sort out some things in it first.

In short, I'm writing an integrator so the eq1 and eq2 parts are contained in a function and hence why I needed the feval.

I suppose keeping eq1 and eq2 seperate is just a personal choice for now to keep things clearer for me.

I don't understand how you can use just one output for fsolve though. I always get the error message saying I too little outputs for my given variables...

You're right about the 2 by 1 vectors being returned although in my original code, k1 (where you have k(1) and k(3)) is a 1 by 2 vector (same with k2) and this will become a 1 by n vector when I sort out my code.

I think I'm gonna leave this fsolve part till I can speak to my lecturer about it... I've already solved this whole thing but I just want to see if I can make it faster...

[Matt Tearle]

As the doc for fsolve (http://www.mathworks.com/help/toolbox/optim/ug/fsolve.html) explains, it needs a function that takes an n-element vector input and returns an n-element vector output, which is why you'll need to join them all into a 4-D system. You could still make eq1 and eq2, then glue them together when defining the function to give to fsolve.

My question about feval was why you're using feval (feval(fun,x,y)) rather than just evaluating the function handle directly (fun(x,y)).

[Brian]

ignore that last deleted comment if you read it, already had out defined as something...