cancel
Showing results for 
Search instead for 
Did you mean: 
cancel
Showing results for 
Search instead for 
Did you mean: 

Community Tip - New to the community? Learn how to post a question and get help from PTC and industry experts! X

Programming of Simpson's rule

AdoS
6-Contributor

Programming of Simpson's rule

Team,

 

Hope you can help me with one more question. 

 

I was solving on-bottom stability and I'm stucked on aproximation of the integral of the function. There are many ways how to solve this Trapezoidal rule, SImpson's rule etc.. so any of those would be great (but not sure how to program any of them). 


Can you please check this spreadsheet and help me. Basicly I'm not sure how to call for even/uneven numbers based on the defined gamma. 


Thank you 

 

AS

1 ACCEPTED SOLUTION

Accepted Solutions

Here for comparison the results with the same value of h using the functional approach.

What you can see is twofold:

1) The approxamation  by simply summing up the rectangles is pretty acurate

2) The 284 special values of omega which you provide have no influence on the outcome - they don't matter.

B2.jpg

View solution in original post

15 REPLIES 15
LucMeekes
23-Emerald III
(To:AdoS)

Why try to write your own numeric integration approximation (Simpson), when you can have Prime do it for you?

Example:

LM_20180720_Integrate.png

Success!
Luc

 

LucMeekes
23-Emerald III
(To:AdoS)

Ah, I see. You don't have functions, you have discrete data/vectors (although I guess that's just your doing, you might as well have omega as a running variable and use the Suu function).

Anyway, in case of vector data, integration becomes summation. See if this helps:

LM_20180720_Integrate1.png

Note! It's not complete, I left it for you to finish. But at least it gives you a clue how to get to even and uneven indexed values...

 

Success!
Luc

Werner_E
24-Ruby V
(To:AdoS)

It wouldn't be hard to implement Simpson but as Luc said its not necessary as Prime would do the job of numeric integration for you if possible.

Unfortunately its not possible with the functions an limits you provide.

You have to evaluate an improper integral (the upper limit is infinity), but the function g(omega), which is part of the integrand, will fail for values of omega greater than  21.8 because of numerical limitation.

Here is a possible workaround using an auxiliary function g_ which is set to zero for omega>21.

B.jpg

B2.jpg

AdoS
6-Contributor
(To:Werner_E)

Thank for the answers guys,

So if I use mathcad prime integration, it would be like this. Looks good, or I did something wrong? :=)Capture232.JPG

Werner_E
24-Ruby V
(To:AdoS)


@AdoS wrote:

Thank for the answers guys,

So if I use mathcad prime integration, it would be like this. Looks good, or I did something wrong? :=)

Not sure. If S.n in your sheet is the same as what is calles S.uu in your template, than it seem so be OK.

You had defines S.uu differently in your sheet and so I used the S.uu you defined.

 

Do you get a result for M(5) or higher values of n?

 

Are you sure that the function which is called S.uu in the picture you show is the same which you call S.n in your sheet and not a function similar to what you call S.uu in your sheet??

AdoS
6-Contributor
(To:Werner_E)

Gents,

 

first sorry for late answer. I had weekend far away from the Internet connection ( it was a bit weird, but relaxing).

 Thank you the answers. Yes I mixed the functions, I think it should be good now, that little trick with w for over 21.8 worked. 

Also for previous solution attached I was not able to calculate M over 5.  

Question is this 21 presents value of w of 21 of position of w 21 which in my case is 0.52?

 

Regards

 

 

Werner_E
24-Ruby V
(To:AdoS)


Question is this 21 presents value of w of 21 of position of w 21 which in my case is 0.52?

 

Sorry, but i don't understand your question.

Values for w which are greater than approx. 21.84 will make your function g(w) fail because the denominator in the expression you calculate exceeds the numerics limit of approx 10^308. The value of g would be lower than 10^-308 and thats the reason its set to zero for w>21 in  the function g_(w)

B0.jpg

But I don't understand what the w in all of your functions and expression should have to do with the vector w you had created. The value with index 21 is 0.52, yes. But the "21" I used as a limit for setting g to 0 is the value of omega itself, not its number in your vector.

BTW - did you typed in that vector omega by hand? All 283 values? As they are equidistant, you could have used something like

B1.jpg

to create that very same vector.

AdoS
6-Contributor
(To:Werner_E)

Hi, 

 

Thank you for the reply, yes I have typed all by hand, thank you for showing me this trick it is fantastic. And thanks for explaining 21.  

 

Acctually I do firmly think I made some error with definition of this w All I wanted is to define values of X axis which shall fill Snn(w) function, and this should be multiplied with G(w) Function and used in integral as is. 

 

That is why I started to write in a vector (first sheet I sent), and use some of approximation rules to calculate M. 

 

jw.JPGCapture1.JPGCapture2.JPG

Werner_E
24-Ruby V
(To:AdoS)

You don't need to define omega to get a plot of your functions (and if, you'd rather define omega as a range -> w:=0.1,0.12..5.74)

Simple use an undefined w in your plot and change the values at the axis as you need.

For each axis you can change the first, second and last value (as highlighted in the plot screenshot below).

Plot looks different as in your picture (hope I got the correct functions to plot), especially the magnitude on the ordinate. You did not use units so we cannot check for consistency here.

B.jpg

AdoS
6-Contributor
(To:Werner_E)

Well idea was not to get plot, rather to get value of M0 which I need to calculate velocity Us (last line in the attached sheet). 

 

But I see the point. Please see attached sheet now with omega defined as a range. nothing works for value of h over 35, but I guess that is due not possibility to iterate. Might assume not so much wave influence at that depth (h is the depth) 

 

Regards

Werner_E
24-Ruby V
(To:AdoS)

For your purpose you don't need to define omega anyway! There is no need to define values. The range you defined is not used (other than in the plot) anywhere in your sheet. You can delete it.

If you need a series of values for a calculation, you should use a vector.

A range should be used just to index the elements of a vector/matrix, for plotting and in a programmed for-loop.

 

For me, the integral already fails for h=33.3 (its the integral, which fails, not any of the solve blocks)

B1.jpg

Maybe it could be an option to do the integration not up to infinity but rather up to a finite value like 10 or 14:

B1.jpg

If thats numerically accaptable, you could use a try and catch construct to fall back to that solution only if the integral up to infinity fails:

B3.jpg

You may apply this principle of course also to the function M(n), but obviously you only need the value of M for n=0.

AdoS
6-Contributor
(To:Werner_E)

Hi, 

 

Thank you for the answer and I'm sorry for so many back and forward. Looks that I poorly define what I want to accomplish. 

 

Yes, I need to use this set of values, they represent frequencies for certain periods of waves. 

 

When I started (somehow got lost in all this), I defined it with vectors, and used vectorization to get function Suu, 

 

Vectorized.JPG

 

But I was not able to solve that M(0) from this function now, and I know it should be simple (not for me). 


Regards

Werner_E
24-Ruby V
(To:AdoS)

As I understand it M(n) is defined by an Integral running from 0 to infinity, not the sum of some discrete values?

So M  has nothing to do with a couple of omega-values from 0.10 to 5.7 but is calculated by spanning over the full continuous range of omega from 0 to infinity.

Of course you could approximate that integral by summing up the product of every value in S.uu with the constant distance of your omega values (0.02). But I see no point in doing so but maybe I am missing something.

B.jpg

Applying a trapezoid rule or Simpson would make the result only very slightly more accurate (it still would only be an integration from 0.12 to 5.74 and not from 0 to infinity). And it still would mean replacing the good approximation from Mathcads solve block functions and its numeric integration by a bold one.

Here for comparison the results with the same value of h using the functional approach.

What you can see is twofold:

1) The approxamation  by simply summing up the rectangles is pretty acurate

2) The 284 special values of omega which you provide have no influence on the outcome - they don't matter.

B2.jpg

AdoS
6-Contributor
(To:Werner_E)

Hi, 

 

Thank you for clarifying this to me sir, you are  a mathcad genius. 

 

I need too look into approach and choose which one suits best, but yes this is what I was looking. 

 

 

 

 

Top Tags