How to calculate steady state with fsolve with white noise in equations ?

I have following equations type

dot(x)=-ax +by^2+n

dot(y)=-cy+dx+w

a,b,c,d are known constants

n and w are gaussian white noise which be defined as n or w =sigma*randn(1,L)

If I dont have n or w .I can easly solve for Steady state using fsolve function or any other function by take dot(x) and dot(y) equal to zero.But how can I handle white noise in steady state solutions?

[John D'Errico]

If a,b,c,d are all known constants, then the steady state solution , thus a limit for x(t) and y(t) as t --> inf. In the absense of noise, you would find this as the solution of the problem

-ax + by^2 == 0

-cy + dx == 0

Now I have no idea if this is your real problem, or just an example one. Usually a real problem is more complicated.

However, as this is a nonlinear problem, then any such nonkinear problem may often have multiple solutions. The one we see here is a simple one. For example, I'll just pick some very simple numbers, and solve the homogeneous problem. Because the problem is so trivial as posed, I can use solve to get all solutions.

xy = solve(-x + y^2 == 0, -y + 2*x == 0)
xy = 
  struct with fields:

    x: [2×1 sym]
    y: [2×1 sym]
    
xy.x
ans =
   0
 1/4
 
xy.y
ans =
   0
 1/2

The nice thing here is solve tells me there are indeed two distinct solutions, whereas fsolve would find one, but not the other, and the solution you find will depend on the starting values.

So there are two solutions to the problem I posed, although without some thought that I don't really want to invest at the moment, I might want to learn if the (0,0) solution is a stable one. My gut tells me it is not a stable one. If so, then any temporary perturbation away from that point, and the solution will drift to the second solution.

Now, you describe a stochastic system, where at any point in time, the differential equation is being hit with a random impulse in some random direction. Is that your question? Are you asking how to analyse such a problem?

If so, then I would suggest that you start with determining if the various potential steady state solutions are stable or not. If one solution is not a stable point, then any perturbation in the system will drive the solution towards the other solution. If you are at a stable solution, then any random impulse will cause the problem to wander around the homogeneously stable solution, in a fuzzy random walk.

For example, suppose you are near/at the (0.25,0.5) homogenous solution for this system? Then imagine your system receives a sudden impulsive perturbation.

xy = solve(-x + y^2 - 0.02 == 0, -y + 2*x + 0.01 == 0)
vpa(xy.x)
ans =
 -0.019194109070750548052986782472964
   0.25919410907075054805298678247296

vpa(xy.y)
ans =
 -0.028388218141501096105973564945927
   0.52838821814150109610597356494593

The steady state points are now subtly different. The system would want to wander in that direction, except that the next time step will again see a random impulse in some random direction.

You might want to characterise the solutions in terms of their stability. You might want to characterize the system in terms of the shape of that fuzzy random cloud you would see around a steady state point.

I would also note that if your problem is less specific, and a,b,c,d are some other values, then the stady state points are:

syms a b c d
xy = solve(-a*x + b*y^2 == 0, -c*y + d*x == 0)
xy = 
  struct with fields:

    x: [2×1 sym]
    y: [2×1 sym]
    
xy.x
ans =
               0
 (a*c^2)/(b*d^2)
>> xy.y
ans =
           0
 (a*c)/(b*d)

But still you need to do more analysis here, beyond this point.

I think you need to define your problem more clearly, since I'm just guessing as to your real question here, postulating what I might want to do, if that were your real question. And since I have no real expertise in differential equations, and certainly no expertise in stochastic problems of this sort, I'll stop here. You probably need to do some reading before you proceed too far, since I'm fairly sure these problems have been studied to a great deal of depth.

Answer 1

[Thiago Henrique Gomes Lobato]

That depends on how you define Steady state. If your derivatives are never actually zero, you will never have a fully perfect steady state, so your problem is impossible if the goal is a perfect steady state. Considering that x and y are spatial coordinates, one way of thinking of this is to find which position will give you the smallest possible variance, which, since your noise is zero mean gaussian, will actually be in the place where:

dot(x)= n

dot(y)= w

or, equivalenty:

-ax +by^2 = 0

-cy+dx+w = 0

That is the same solution you would become had you solved the problem completely ignoring the noise. This is related to the Maximum likelihood estimation and it basically says "Considering that my data has zero mean Gaussian noise, which is the most likely solution between all solutions?", and this is the one where your mean "error" is zero and all that's left are the noise fluctuations.

[DEVENDER garg]

Thanks for your reply!.I think it will give correct solution if we look at value of x and y after simulation.But I think if we look at value of x^2 or y^2 after simulation it can heavly depend on w*w or n*n. one way,I am thinking of using for loop each one for particular value of vector randn(1,L).Then average the the output result.But in principal it does not look correct to me.

Another one is using differntial equation direclty with time and noise.And manually looking for steady state which has error in the range of gaussian noise as you suggested.