Advanced Search

hajnayeb

15 Reputation

3 Badges

6 years, 80 days
Contact hajnayeb

MaplePrimes Activity


These are replies submitted by hajnayeb

Re:...

hajnayeb 15
August 10 2020
0 0

@mmcdara 

Hi! I changed your code in a way that I can use it for finding the variance of the response of an SDOF system to a random noise.

I like to have a shite noise with a variance of 4 and zero mean. It seems I don't know how to set these numerical values.

Plus, the code seems simple but take a long time to run.

 

SDOF-Linear.mw
 

> 

restart; with(LinearAlgebra); with(plots); with(Statistics)
> 

KG := proc (X, Y, psi, sigma) local NX, DX, K, mu, k, y; NX := numelems(X); DX := `<,>`(seq(Vector[row](NX, X), k = 1 .. NX))^2; K := proc (sigma, psi) options operator, arrow; evalf(`~`[`*`](sigma^2, `~`[exp](-`~`[`^`](`~`[`/`](DX-LinearAlgebra:-Transpose(DX), psi), 2)))) end proc; mu := add(Y)/NX; k := proc (t, sigma, psi) options operator, arrow; evalf(convert(`~`[`*`](sigma^2, `~`[exp](-`~`[`^`](`~`[`/`](`~`[`-`](t, X), psi), 2))), Vector[row])) end proc; y := mu+k(t, sigma, psi).(1/K(sigma, psi)).convert(`~`[`-`](Y, mu), Vector[column]); return y end proc
> 

t_end := 2; NX := 500; d_t := evalf(t_end/NX)

500

  

0.4000000000e-2

 (1)
> 

T := [seq(0 .. t_end-d_t, d_t)]
> 

Nb_traj := 10
> 

z := .2; wn := 10; m := 1

.2

  

10

  

1

 (2)
> 

k := wn^2*m; c := 2*z*sqrt(k*m); mxn := Array(1 .. Nb_traj, 1 .. 2)

100

  

4.0

 (3)
> 

xIC := 0; vIC := 0

0

  

0

 (4)
> 

for traj to Nb_traj do F := convert(Statistics:-Sample(Normal(0, 4), NX), list); X := dsolve({m*(diff(x(t), `$`(t, 2)))+c*(diff(x(t), t))+k*x(t) = KG(T, F, .1, 1), x(0) = xIC, (D(x))(0) = vIC}, numeric); pl || traj := odeplot(X, [t, x(t)], t = 0 .. 10, color = ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])) end do
> 

display(seq(pl || traj, traj = 1 .. Nb_traj))

 
> 

``
> 

for traj to Nb_traj do F := convert(Statistics:-Sample(Normal(0, 4), NX), list); X := dsolve({m*(diff(x(t), `$`(t, 2)))+c*(diff(x(t), t))+k*x(t) = KG(T, F, .1, 1), x(0) = xIC, (D(x))(0) = vIC}, numeric); dispx := proc (u) options operator, arrow; subs(X(u), x(t)) end proc; velx := proc (u) options operator, arrow; subs(X(u), diff(x(t), t)) end proc; xn := [seq(dispx(i), i = 0 .. t_end, 2)]; mxn[traj, 1] := Mean(xn); mxn[traj, 2] := StandardDeviation(xn)^2 end do
> 

``
> 

mxn

Array(%id = 18446746202355763790)

 (5)
> 

Mean(mxn)

Vector[row](%id = 18446746202355810902)

 (6)
> 

``
> 

``
> 

``


 

Download SDOF-Linear.mw

 

Monte-Carlo...

hajnayeb 15
August 06 2020
0 0

@mmcdara 

Yes, vallidating the exact results seems strange but in articles, usually the reviewers ask for validation to chech the correctness of the numerical results. 

I would liike to assess the mean and variance of X(t) for any Gaussian white-noise input with MCM?

Using what you sent, I can get several outputs and find the averaged values. I think I have all the necessary material. :)

Frequency-domain method...

hajnayeb 15
August 05 2020
0 0

@mmcdara 

I know that an easy way for linear systems is using the frequency response function, but I am looking for a validation technique (Monte-Carlo simulation).

Thanks...

hajnayeb 15
August 05 2020
0 0

@mmcdara 

It seems that it is the algorithm I was looking for.

I have to implement it myself. I will get back to you as soon as I am able to derive everything myself.

BTW, attached is a simulink file, which simulates the same system, but it does not give accurate results in case of MDOF systems. 

Thanks a lot!

 

Re: What do you mean by "f(t) is a rando...

hajnayeb 15
August 05 2020
0 0

@mmcdara 

Thanks for your reply!

As you mentioned, I am interested to apply a Gaussian White Noise, with a specific mean and variance, to a mehanical system and determine the mean and variance of the response. In simulink, there is a "Band-limited white noise" block for that. We can use mean and var commands to determine the output mean and variance, too.

Because of some problems of Matlab, I would like to implement that in Maple.

Thanks again!

The excitation you created is not random...

hajnayeb 15
August 05 2020
0 0

@nm 

Thanks for your answer.

I checked the code. The excitation you created as F(t) is exponential not random. If you check the following version of your code, even with zero ICs, it has nonzero periodic response. I plotted F(t) too.

 

Random_testfile.mw
 

> 

restart
> 

f := Statistics:-RandomVariable(Normal(0, .1)):
> 

plot(Statistics:-PDF(f, t), t = -5 .. 5);

 
> 

ode := m*(diff(x(t), `$`(t, 2)))+c*(diff(x(t), t))+k*x(t) = F(t);

m*(diff(diff(x(t), t), t))+c*(diff(x(t), t))+k*x(t) = F(t)

 (1)
> 

ic := x(0) = 0, (D(x))(0) = 0;

x(0) = 0, (D(x))(0) = 0

 (2)
> 

m := 1;

1

  

.3

  

10

 (3)
> 

F := proc (t) options operator, arrow; Statistics:-PDF(f, t) end proc;

proc (t) options operator, arrow; Statistics:-PDF(f, t) end proc

 (4)
> 

sol := dsolve([ode, ic], x(t));

x(t) = -(5/3991)*exp(-(3/20)*t)*sin((1/20)*3991^(1/2)*t)*exp(-1991/40000)*(exp(-((3/40000)*I)*3991^(1/2))*3991^(1/2)*erf((1/400)*2^(1/2)*(I*3991^(1/2)-3))*Pi^(1/2)-(200*I)*exp(-((3/40000)*I)*3991^(1/2))*exp(-(1/80000)*(I*3991^(1/2)-3)^2)*2^(1/2)-exp(((3/40000)*I)*3991^(1/2))*erf((1/400)*2^(1/2)*(I*3991^(1/2)+3))*3991^(1/2)*Pi^(1/2)+(200*I)*exp(((3/40000)*I)*3991^(1/2))*exp(-(1/80000)*(I*3991^(1/2)+3)^2)*2^(1/2))/Pi^(1/2)+exp(-(3/20)*t)*cos((1/20)*3991^(1/2)*t)*(((5/3991)*I)*exp(-((3/40000)*I)*3991^(1/2))*erf((1/400)*2^(1/2)*(I*3991^(1/2)-3))*exp(-1991/40000)*3991^(1/2)+((5/3991)*I)*exp(((3/40000)*I)*3991^(1/2))*erf((1/400)*2^(1/2)*(I*3991^(1/2)+3))*exp(-1991/40000)*3991^(1/2))-(5/3991)*exp(-(3/20)*t-1991/40000)*3991^(1/2)*(exp(-((3/40000)*I)*3991^(1/2))*(I*cos((1/20)*3991^(1/2)*t)-sin((1/20)*3991^(1/2)*t))*erf((1/400)*2^(1/2)*(I*3991^(1/2)+2000*t-3))+exp(((3/40000)*I)*3991^(1/2))*erf((1/400)*2^(1/2)*(I*3991^(1/2)-2000*t+3))*(I*cos((1/20)*3991^(1/2)*t)+sin((1/20)*3991^(1/2)*t)))

 (5)
> 

plot(rhs(sol), t = 0 .. 20)

 
> 

plot(F(t), t = 0 .. 4)

 
> 

``


 

Download Random_testfile.mw

Thanks...

hajnayeb 15
October 03 2019
0 0

@Kitonum 

I changed it as you wrote, now it works.

Page 1 of 1