Fsolve
in [Matlab]
|
From: Walter Roberson on 29 Mar 2010 19:51 rabbahs wrote: > Walter Roberson <roberson(a)hushmail.com> wrote in message > <hoinjr$r39$1(a)canopus.cc.umanitoba.ca>... >> rabbahs wrote: >> It means what I said before: that I would not recommend trying to >> solve the problem with a minimizer. I showed the form of the exact >> solution of the problem, and if you have an exact solution, you do not >> need a minimizer: you just plug in the numbers, >> I could recompute the exact solution (sorry, I deleted it yesterday... >> on the other hand I might well be able to recover it from my Time >> Machine backups) and email it to you in Maple format, but it was about >> 12 megabytes. Actually it was 19 megabytes of incomprehensible masses of symbols. > If u recover the > solution then kindly let me know I will not send the final solution because it is just too big. Instead, I will show how to get the solution in Maple. I highly recommend using the command-line version of Maple for this purpose. First, reduce down to the 3 equations in a, b, d, like we discussed before. I will refer to this as eqmat3. eqmat3sol := solve(eqmat3,[a,b,d]) assuming real; OR// eqmat_sol := RealDomain:-solve(eqmat3,[a,b,d]); If you are using the command line version of Maple, it will print out a single solution of the form [[a = mess#1 involving %1, b = %1, d = mess#2 involving %1]] and then it will print out %1 := mess#3 This means that Maple has identified a significant common sub-expression to several terms, and has broken it out so as to make the output easier to read. If you were using one of the GUI versions of Maple, you would have to go through a bit of trouble to get Maple to print things out with this common sub-expression made clear. The common sub-expression will be a RootOf() call close to 60 lines long that sets out a quartic equation. There will be four roots to the quartic, and they will all be terribly long and messy; you will have a fair bit of trouble understanding the solution if you proceed to ask Maple to find all four roots of the quartic. And it will take a long time, too. So we will step around the issue a bit to make things faster. Now, command: P1val := %1: This isolates the RootOf() into a variable. Now construct a variant of the solution matrix rewritten with a symbolic variable for this RootOf: eqmat3_solP1 := subs(P1val=P1,eqmat3_sol): Once we have found candidate solutions for the RootOf(), we will go back and substitute them in individually for P1, and ensure that a and d come out non-negative (b is the same as the root, and since b must non-negative as well, we will only select non-negative roots.) If a and d pass, then we go back and substitute the a, b, and d into the other 5 variables that have been rewritten in terms of a, b, and d, as shown earlier in the discussion. The next step towards finding the solution is formally necessary according to the Maple documentation, but I've never seen a problem if it isn't used. It certainly doesn't hurt and probably makes things faster, so we might as well do it anyhow: Roc := collect(op(P1val),_Z): C4v := factor(coeff(Roc,_Z,4)): C3v := factor(coeff(Roc,_Z,3)): C2v := factor(coeff(Roc,_Z,2)): C1v := factor(coeff(Roc,_Z,1)): C0v := factor(coeff(Roc,_Z,0)): This isolates the coefficients of the _Z terms of the RootOf. Now the time saving step: instead of finding the RootOf the big messy equation, we will find the RootOf a generic quartic, and substitute in the above coefficients afterwards. G4 := RootOf(C4*_Z^4+C3*_Z^3+C2*_Z^2+C1*_Z+C0): G4S := [allvalues(G4)]: I pre-computed G4S and saved it away so that I don't have find the solutions over and over. tmpproc := unapply(G4S[1],[C0,C1,C2,C3,C4]): G4S1proc := codegen[optimize](tmpproc): tmpproc := unapply(G4S[2],[C0,C1,C2,C3,C4]): G4S2proc := codegen[optimize](tmpproc): tmpproc := unapply(G4S[3],[C0,C1,C2,C3,C4]): G4S3proc := codegen[optimize](tmpproc): tmpproc := unapply(G4S[4],[C0,C1,C2,C3,C4]): G4S4proc := codegen[optimize](tmpproc): After this set of Maple commands, G4S*proc are procedures for efficiently computing the individual quartic roots given the coefficients. Again, you may wish to save these procedures away for future use. Now: P1sol1 := G4S1proc(C0,C1v,C2v,C3v,C4v): P1sol2 := G4S2proc(C0,C1v,C2v,C3v,C4v): P1sol3 := G4S3proc(C0,C1v,C2v,C3v,C4v): P1sol4 := G4S4proc(C0,C1v,C2v,C3v,C4v): These are now the individual roots for P1, the RootOf(), which is also the value for b. If you want to proceed symbolically from here, you could construct entire solution sets using (e.g.) eval(eqmat3_solP1, P=P1sol1); but you are going to get messy. I would recommend this as being about the point that you switch from symbolic to numeric. Code as far as the P1sol* in terms of the constants c1(1), c1(2) and so on (which I refer to as c11, c12 and so on here), and then code the eqmat3_solP1 in terms of the constants and P1, and then code the rest of the variables in terms of the a, b, and d that you get out. Only bother to plop in the values of P1 that came out nonnegative and real; only bother to continue on to the extended equations if a and d also come out nonnegative and real; and, of course, afterwards sanity-check everything to be sure that all the variables came out non-negative and real. With this construction, if any of the variables come out negative or complex, then there is no solution for that set of c1 constants. If you are going to be using a number of different c1 constant sets, then you might want to take the extra step of: tmpproc := unapply(P1sol1,[c11, c12, c13, c14, c15, c16, c17, c18]): P1sol1proc := codegen[optimize](tmpproc): tmpproc := unapply(P1sol2,[c11, c12, c13, c14, c15, c16, c17, c18]): P1sol2proc := codegen[optimize](tmpproc): tmpproc := unapply(P1sol3,[c11, c12, c13, c14, c15, c16, c17, c18]): P1sol3proc := codegen[optimize](tmpproc): tmpproc := unapply(P1sol4,[c11, c12, c13, c14, c15, c16, c17, c18]): P1sol4proc := codegen[optimize](tmpproc): (The roots should come out pretty similar, actually -- just like a quadratic has two solutions because it has a basic formula with one place to choose +/-, a quartic has four solutions because it has a basic formula with two places to choose +/- .) |