Minimizing a mutivariable symbolic function

My code reads as follows:

syms u P Q
U = sym('u', [101 1]);
P(1,1) = 1;
Q(1,1) = 0; 
for k = 1:100
    P(k+1,1) = .05*(-.3*P(k,1) +.65*Q(k,1)+U(k,1))+P(k,1);
    Q(k+1,1) = .05*(-.65*Q(k,1) + .3*P(k,1)-U(k,1))+Q(k,1);
end
A = P-.85;
B = Q-.15;
Z = [A; B];
J = norm(Z)^2 + norm(U)^2;

I simply want to minimize the function J, which is a function of the symbolic variables u1,...,u100(elements of the vector U). However, for some reason I am not able to implement the script as a function (possibly because of the symbolic terms) and so I simply have the file saved as "Uvector". But when I attempt to run " a = fminsearch(@Uvector, zeros(101,1))" to find the minimizer, an error reads "Attempt to execute SCRIPT Uvector as a function". Any suggestions on how to save the file so that J is executed as a function or another command that may minimize the function?

[Walter Roberson]

To check: your desired output would be a function of 100 variables, u1 to u100, that you could then put through something like fmincon() ?

Answer 1

[Walter Roberson]

Do you have access to the Symbolic Toolbox? If so then you might want to take advantage of the recurrence equation solver built into MuPAD; see http://www.mathworks.com/help/symbolic/mupad_ref/rec.html

I find that

P(n) = sum(20^(-1-2*n+2*n0)*381^(n-n0) * U(n0-1), n0 = 2 .. n) + 13/19 + (800/2413) * (381/400)^(-n)
Q(n) = -(sum(20^(-1-2*n+2*n0)*381^(n-n0) * U(n0-1), n0 = 2 .. n)) + 6/19 - (800/2413) * (381/400)^n

(Note: This is not normal MATLAB notation; sum() is part of the Symbolic Toolbox)

Notice that each p(n) and q(n) involve all of the U(1:n). Because of this, the 2-Norm expands fairly quickly when calculated symbolically. You are probably not going to have enough time and memory to calculate the symbolic J function out to 100.

I suspect you are going to find it a lot more efficient to make the calculations numerically for each set of U values. Calculating symbolically even just to 5 takes a fair bit of time and memory.

=====

As I have already worked it out, I include the Maple version of the code, below. N would be 100 in your original problem, but you are not going to get even a fraction of that high without running out of memory:

N := 5:
with(LinearAlgebra):
P := Vector(N+1):
Q := Vector(N+1):
U := Vector(N, symbol = 'u'):
P[1] := 1:
Q[1] := 0:
for k to N do
  P[k+1] := 1/20*(-(3/10)*P[k]+(65*(1/100))*Q[k]+U[k])+P[k]:
  Q[k+1] := 1/20*(-(65*(1/100))*Q[k]+(3/10)*P[k]-U[k])+Q[k]:
end do:
A := `~`[`-`](P, 85*(1/100)):
B := `~`[`-`](Q, 15*(1/100)):
Z := `<,>`(A, B):
J := VectorNorm(Z, 2, conjugate = false)^2 + VectorNorm(U, 2, conjugate = false)^2;

And the recurrence, again in Maple syntax:

rsolve({p(1) = 1, p(K+1) = 1/20*(-(3/10)*p(K)+(65*(1/100))*q(K)+u(K))+p(K), q(1) = 0, q(K+1) = 1/20*(-(65*(1/100))*q(K)+(3/10)*p(K)-u(K))+q(K)}, {p(n),q(n)})

You can build the J function iteratively: the n'th additional term adds p(n)^2 + q(n)^2 + u(n)^2 to the running total.