fsolve, time integration,

Hallo everybody,

I'm new here and i need help using fsolve! I would like to solve the following equation: All values are known except x

 F = @(x) M*x + D*UTDOT + K*UT - Fext;  

 a = fsolve(F,0);

where

ALGAT = ((1-af)/(1-am))*x + (af/(1-am))*a - (am/(1-am))*alga;

UT = u + h*v + ((h^2)/2)*( (1-2*beta)*alga + 2*beta*ALGAT );

UTDOT = v +h*( (1-gamma)*alga + gamma*(ALGAT) );

if i write my problem the following way, i'm able to obtain a solution

    F = @(x) M*x+D*(v +h*( (1-gamma)*alga + gamma*(((1-af)/(1-am))*x + ...
        (af/(1-am))*a - (am/(1-am))*alga) ))...
        +K*(u + h*v + ((h^2)/2)*( (1-2*beta)*alga + ...
        2*beta*((1-af)/(1-am))*x + (af/(1-am))*a - (am/(1-am))*alga ))-Fext;
    a = fsolve(F,0);

But i would like to solve the system in this form

 F = @(x) M*x + D*UTDOT + K*UT - Fext;      %M,D,K,Fext known values
 a = fsolve(F,0);

How can handle UTDOT and UT to F ???

Thank you for your help!

Answer 1

[Star Strider]

I cannot run your code, but if ALGAT is a funciton of ‘x’, it and every term that uses it also need to be.

See if this works:

ALGAT = @(x) ((1-af)/(1-am)).*x + (af/(1-am))*a - (am/(1-am))*alga;

UT = @(x) u + h*v + ((h^2)/2)*( (1-2*beta)*alga + 2*beta.*ALGAT(x) );

UTDOT = @(x) v +h*( (1-gamma)*alga + gamma.*(ALGAT(x)) );

F = @(x) M*x + D.*UTDOT(x) + K.*UT(x) - Fext;      %M,D,K,Fext known values

[Stefan Holzinger]

Hallo Star Strider,

Thank you for your help. Your solution works! :) But if i run the whole code:

 t0 = 0;
 tend = 5;
 h = 0.01;

 M = 0.2;
 D = 0.5;
 K = 15;
 Fext = 0;

 p0 = 0.2;
 v0 = 0;
 a0 = 0;

 rinf  = 0.98;  
 am    = (2*rinf-1)/(rinf+1);
 af    = rinf/(rinf+1);
 beta = (1/4)*(1-am+af)^2;
 gamma = (1/2)-am+af;

 t = zeros((tend-t0)/h,1);

 for i = 1:length(t)

    t(1) = t0;
    u(1) = p0;
    v(1) = v0;
    a(1) = 0;
    alga(1) = 0;

    [M,D,K,Fext,dCdu,C] = F_MBS_TRY1();

    ALGAT = @(x) ((1-af)/(1-am)).*x + (af/(1-am))*a(i) - (am/(1-am))*alga;
    UT    = @(x) u(i) + h*v(i) + ((h^2)/2)*( (1-2*beta)*alga + 2*beta.*ALGAT(x) );
    UTDOT = @(x) v(i) +h*( (1-gamma)*alga + gamma.*(ALGAT(x)) );

    F = @(x) M*x + D.*UTDOT(x) + K.*UT(x) - Fext;
    a(i+1) = fsolve(F,0);

    alga(i+1) = ((1-af)/(1-am))*a(i+1) + (af/(1-am))*a(i) - (am/(1-am))*alga(i);    
    u(i+1)    = u(i)+h*v(i)+((h^2)/2)*((1-2*beta)*alga(i)+2*beta*alga(i+1));
    v(i+1) = v(i)+h*((1-gamma)*alga(i)+gamma*alga(i+1));
    t(i+1) = t(i)+h;

 end

fsolve returns the following Warning:

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead.

I checked a few sites concerning this topic but i still don't know what to do with this error. Do have an idea how i can get rid of this warning?

[Star Strider]

My pleasure!

There is nothing wrong. It is a ‘warning’, not an ‘error’. The fsolve function is telling you that it is switching to the Levenberg-Marquardt algorithm. It should not affect the result of the optimisation.

[Stefan Holzinger]

Hallo Star Strider,

Ah I see, OK thanks! It seems to me that you are Matlab Pro :) There is one last thing which causes me some headache. I wrote in the above mentioned time integration context my own newton-raphson code. If i run the code that i posted earlier, matlab returns the following error message:

 Assignment has more non-singleton rhs dimensions than non-singleton subscripts

 Error in newton (line 33)
        DF(:,i) = ( F(x0+h*E(:,i)) - F(x0-h*E(:,i)) )/(2*h);

 Error in alpha1 (line 31)
    a(i+1) = newton(F,0);

This is my newton code:

 function [x,f,loopcount,ConditionNumbDF] = newton(F,x0)

 error = 1;
 tolerance = 1.e-8;  

 h0 = (eps)^(1/3);
 n = length(x0);
 DF = zeros(n);
 E = eye(n);

 loopcount = 0;
 while error > tolerance

    % Jacobian 
    for i = 1:n
        h = max(1,abs(x0(i)))*h0;
        DF(:,i) = ( F(x0+h*E(:,i)) - F(x0-h*E(:,i)) )/(2*h);
    end

    x = x0 - DF\F(x0)';
    error = norm(x-x0);
    x0 = x;

    loopcount = loopcount + 1;
    f(:,loopcount) = F(x0);
    ConditionNumbDF(loopcount) = cond(DF);
end
end

Somehow i can't find the mistake. Do you have any ideas?

[Star Strider]

I can’t run your code, so executing only that loop:

F = @(x) x.^2 - 1;

x0 = 0;
h0 = (eps)^(1/3);
n = length(x0);
DF = zeros(n);
E = eye(n);

for i = 1:n
    h = max(1,abs(x0(i)))*h0;
    DF(:,i) = ( F(x0+h*E(:,i)) - F(x0-h*E(:,i)) )/(2*h);
end

it runs for me without error. If ‘x0’ has more than one element, it must be a column vector for your code to work, so adding this line early in your code will eliminate that problem:

x0 = x0(:);

That will create a column vector out of any vector of ‘x0’. That is not a problem here with only one parameter.

[Stefan Holzinger]

Hey Star Strider, thank you for your help. Your solution works. :)

In the meanwhile another problem appeared. I tried to solve it but until now i failed. Maybe you have an idea how to do it. Heres the Problem.

I would like to solve the equation F for x :

 ALGAT = @(x) ((1-af)/(1-am)).*x + (af/(1-am))*a - (am/(1-am))*alga;
 UT    = @(x) u + h*v + ((h^2)/2)*( (1-2*beta)*alga + 2*beta.*ALGAT(x) );
 UTDOT = @(x) v +h*( (1-gamma)*alga + gamma.*(ALGAT(x)) );

 [M,D,K] = fun(UT(x),UTDOT(x),t);

 F = @(x) M(UT(x))*x + D(UT(x),UTDOT(x))*UTDOT(x) + K*UT(x) = 0

 x0 = zeros(...)
 solu = newton(F,x0)

The Problem is now: M,D,K are defined within the function fun and they are depend on UT and UTDOT.

My question is: How can i tell matlab within F that M,D and K are dependent of UT(x) and UTDOT(x)??

Answer 2

[Stefan Holzinger]

Hallo Star Strider,

could you maybe take a look a this Problem of mine? Of interest would be my last post in this question. Thank you for your help!

https://de.mathworks.com/matlabcentral/answers/298154-functions-retuning-functions-fsolve

The key words to sind the side would be: functions retuning functions, fsolve