Solving Trigonometric Function Equations with Time varying Parameters by MATLAB

Hello, my goal is to solve the equations about cos (x) and sin (x), plot the solution results into a curve, and the Lac in the equation will change with time.

I tried to use fsolve to solve the equation, but the result is only a fixed value, not a value that changes with time in theory. Here is my code. I look forward to receiving your guidance

%The function I created is as follows

function q=myfun(p)

x=p;

t=0:0.001:5;

Lac=0.91+0.01*sin(t);

q=2 * ( 0.0822*cos(x) - 0.2838*sin(x)) + Lac.^2 - 0.9209;

end

%Here is the command that I call the function

[x]=fsolve(@myfun,[0])

Answer 1

[Davide Masiello]

You were almost there

t = (0:0.001:5)';
x = fsolve(@(x)myfun(x,t),zeros(size(t)));

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

plot(t,x)
xlabel('t')
ylabel('x')

function q = myfun(x,t)
Lac = 0.91+0.01*sin(t);
q   = 2*(0.0822*cos(x)-0.2838*sin(x))+Lac.^2-0.9209;
end

[LIULIYAN]

Thank you very much for your reply. Your answer is very helpful to me. May I ask you one more question.

Lac is not 0.91+0.01 * sin (t), but a variable solved by SIMULINK model (see figure). My ultimate goal is to calculate radian x in real time in SIMULINK according to the real-time changes of Lac.

I tried to use the code you provided to establish the S-function function ("Positivesolution" in the figure). The code is as follows. I found that the running speed was very slow. I ran SIMULINK for 30S (Fixed step size=0.001), but it took several decades to complete the execution, which obviously did not meet the real-time requirements.

Please ask whether the following code can be optimized to improve the real-time performance of SIMULINK？

Or, according to the method in the initial question, how to use the data in simout1.data to assign values to Lac to calculate radian x offline？

Look forward to your guidance, thank you again

function[sys,x0,str,ts] = Positivesolution(t,x,u,flag)
switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes;
    case 1,
        sys=mdlDerivatives(t,x,u);
    case 3,
        sys=mdlOutputs(t,x,u);
    case {2,4,9}
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function[sys,x0,str,ts]=mdlInitializeSizes
global cij bj c
sizes = simsizes;
sizes.NumContStates=1;
sizes.NumDiscStates=0;
sizes.NumOutputs=1;
sizes.NumInputs=1;
sizes.DirFeedthrough=1;
sizes.NumSampleTimes=0;
sys=simsizes(sizes);
x0  = [0.873];
str = [];
ts  = [];
function sys=mdlDerivatives(t,x,u)
sys(1)=x(1);
function sys=mdlOutputs(t,x,u)
Lac=u(1);
%q=2 * ( 0.0822*cos(x) - 0.2838*sin(x)) + Lac.^2 - 0.9209;
x = fsolve(@(x)2 * ( 0.0822*cos(x) - 0.2838*sin(x)) + Lac.^2 - 0.9209,zeros);
sys(1)=x(1);

[Davide Masiello]

If I understand correctly, you want to change the code so that Lac is not a defined function but an array of values calculated in some other code block before. Am I right?

[LIULIYAN]

Yes, sir

[Davide Masiello]

Well then, instead of passing t to the function and defining Lac(t), you just pass it the array of values like below

Lac = 0.91+0.01*sin((0:0.001:5)') % Lac is now an array of doubles

Lac = 5001×1

0.9100 0.9100 0.9100 0.9100 0.9100 0.9100 0.9101 0.9101 0.9101 0.9101

x = fsolve(@(x)myfun(x,Lac),zeros(size(Lac)));

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

plot((0:0.001:5)',x)
xlabel('t')
ylabel('x')

function q = myfun(x,Lac)
q   = 2*(0.0822*cos(x)-0.2838*sin(x))+Lac.^2-0.9209;
end

[LIULIYAN]

It's my pleasure to communicate with you. Your answer has solved my doubts. Thank you very much.

I will continue to study how to solve this equation in real time under the premise of ensuring the running speed of the program in SIMULINK

Answer 2

[John D'Errico]

Way better than using fsolve is to compute the analytical solution. That is, if t is the time varying parameter, then we have

syms t x

Lac=0.91+0.01*sin(t);
q = 2 * ( 0.0822*cos(x) - 0.2838*sin(x)) + Lac.^2 - 0.9209;

Now, we wish to solve q==0, for x, but as a function of t.

xsol = solve(q == 0,x,'returnconditions',true)

xsol = struct with fields:
             x: [2×1 sym]
    parameters: k
    conditions: [2×1 sym]

So there are infinitely many solutions, valu=id for any integer value of k.

xsol.x

ans = 

That may look a bit messy, but we can first take only the primary solution, where k == 0. Agaoin, add any integer multiple of 2*pi to these solutions.

syms k
xprimary = subs(xsol.x,k,0)

xprimary = 

Note there appear to be complex terms in this. We might want to plot the solutions, looking to see if and where imaginary terms enter in. When we do so, we will see the imaginary terms are essentially always zero. So they ended up canceling out always.

fplot(imag(xprimary(1)),[0,5],'b')

fplot(real(xprimary(1)),[0,5],'b')
hold on
fplot(real(xprimary(2)),[0,5],'r')
ylim([-3,1])
grid on
xlabel t
ylabel x(t)

These solutions are in terms of radians. If you want degrees, multiply by 180/pi. but since you used sin and cos in your equaritnis, I can only assume you want to see it in terms of radians.

Finally, IF you want a functional form you can call, then just use matlabFunction. Thus...

x1_t = matlabFunction(real(xprimary(1)));
x2_t = matlabFunction(real(xprimary(2)));

The positive solution was the second one. As you see, we can evaluate it directly for ANY value of t.

x2_t(3)

ans = 0.1286

Thus 0.1286 radians.

[LIULIYAN]

Thank you. Your answer is very helpful to me.