numerical integration and solving for limit

I have the following equation:

g(x) I have as a matrix in an array. It has an x and y column. So I want to take g(x) at a given poin, multiply by 1/x, then integrate, and find where the above statement is true. How can I write a script to solve for X_M for the above integral to be true?

Answer 1

[Star Strider]

You have told us nothing about ‘g(x)’. Assuming the integral of ‘g(x)’ is monotonically increasing at least until its integral is equal to 1 (if it’s periodic or has other pecularities, you will need another approach), try this:

x = linspace(1, 10);                            % ‘x’
y = rand(size(x));                              % ‘g(x)’

vint = pi*sqrt(2)/3 * cumtrapz(x, y./x);
X_M = interp1(vint, x, 1);

Experiment to get the result you want.

[Benjamin]

g(x) is not monotonically increasing. It fluctuates up and down. My curent code looks like:

 [~,sheet_name]=xlsfinfo('C:filename.xlsx');
 for k=1:numel(sheet_name)
   data{k}=xlsread('C:filename.xlsx',sheet_name{k});
   % end
for i=1:1
    plot(data{1, i}(:,2),data{1, i}(:,4));
    C = data{1, i}(:,4)./data{1, i}(:,2);
    grid on;
    hold on;
end

The line I most recently added is C which divides g(x) by x. But when I change this to C(i) and change the loop to 1:30, I get that the left and right sides have different number of elements.

Sorry if this is confusing. Perhaps I could just create an array with the integration and then I can visually see where the solution equals 1? I've attached an example of g(x).

Note that here:

data{1, i}(:,2) is x
data{1, i}(:,4) is y

How can I do this in a loop?

I have something like this below, but how to do this for all in the loop?

Below seems to work. How can I store in array and loop through all of them?

    vint = pi*sqrt(2)/3 * cumtrapz(data{1, i}(:,2), data{1, i}(:,4)./data{1, i}(:,2));
    X_M = interp1(vint, data{1, i}(:,2), 1);

    

[Star Strider]

I would be less confused if I had ‘data’ to work with.

My code divides ‘g(x)’ (whatever it is) by‘ ‘x’ using element-wise operations in the cumtrapz call. A separate element-wise division before using my code is not necessary.

Try my code with ‘data’ or post ‘data’ here so I can see it and work with it. My posted code is the best I can do without it.

[Benjamin]

here is data, can you see it? I think your code actually works for what I want. Look at my last 2 lines in my reply to your comment. How do I extend that to allow for the loop?

[Benjamin]

I think your code works, I just need it to work inside my loop. So then if my loop goes from 1:30, I will have a colmn vector of 30 xM values. Note that the length of the vector may be different for each loop. Currently it just stores the last value of XM computed in the loop

[Star Strider]

I was able to approximate your function with this code, and got a finite result for :

x = linspace(1, 7);                            % ‘x’
y = -25*cos(4*x) .* exp(-3*x) + 1;             % ‘g(x)’

vint = pi*sqrt(2)/3 * cumtrapz(x, y./x);
X_M = interp1(vint, x, 1);

producing:

X_M =
     2.0124

You should get something similar, because your ‘g(x)’ function converges to 1.

[Benjamin]

I think your other code worked just fine. But if I change my loop size, it only records the last value. How do I record all values so they dont get overwritten with each loop?

[Benjamin]

I have this below, but X_M only records the last X_M in the loop (iteration 17). How can I store all of them?

for i=1:17
    plot(data{1, i}(:,2),data{1, i}(:,4));
    C = data{1, i}(:,4)./data{1, i}(:,2);
    
    vint = pi*sqrt(2)/3 * cumtrapz(data{1, i}(:,2), data{1, i}(:,4)./data{1, i}(:,2));
    X_M = interp1(vint, data{1, i}(:,2), 1);
    grid on;
    hold on;
end

[Star Strider]

I’ve not tried your loop because I have to go to Internet Explorer to download .mat files (Firefox has serious problem with them, and it’s even difficult to find where Firefox puts them).

However looking at your loop code, it appears all you need to do is to subscript ‘X_M’:

X_M(i) = interp1(vint, data{1, i}(:,2), 1);

It’s a scalar quantity, so that should work.

[Benjamin]

oh yea, that's it! thanks! so the interp1 just finds where the integral is equal to 1?

[Star Strider]

As always, my pleasure!

Yes.

The assumption is that the integral (produced here by cumtrapz) is monotonically increasing. (If it isn’t, a solution is still theoretically possible, although the code is more complicated.)

[Benjamin]

Thanks again for the help! I had a quick question, and I can make this a new question if I need to, but it seems related -

How can I plot the value of X_M for each plot? Each X_M should have a y associated with it. For each line plotted on the graph, I want to plot this data point with a symbol. This will help me better visualize where X_M falls on each plot. Is this doable?

Basically, wherever the X value of XM is for each line in the loop, I want to plot a symbol at the value of XM on that line, if that makes sense.

[Star Strider]

For the cell in the code you posted, try this:

D = load('data.mat');
data = D.data;

figure
hold all
for i=1:17
    plot(data{1, i}(:,2),data{1, i}(:,4));
    C = data{1, i}(:,4)./data{1, i}(:,2);
    vint = pi*sqrt(2)/3 * cumtrapz(data{1, i}(:,2), data{1, i}(:,4)./data{1, i}(:,2));
    X_M(i) = interp1(vint, data{1, i}(:,2), 1);
    Y_M(i) = interp1(data{1, i}(:,2), data{1, i}(:,4), X_M(i));
    plot(X_M(i), Y_M(i), '+k')
    grid
end
hold off
xlim([1 5])

This uses essentially the same idea (an interp1 call) to calculate the y-value (that I chose to call ‘Y_M’ for consistency) associated with it, and plots it as a black ‘+’. This also saves the values of ‘Y_M’ as a vector for subsequent use. The xlim call eliminates the flat part of the plot in order to make thed plotted ’+’ markers more visible. (I apologise for the delay — I had to go back and remember what we were doing.)

[Benjamin]

wow, thanks! works perfectly! I'll have to dig into this to see what you did. I may ask a question later if I can't figure it out. you have been super helpful!

[Star Strider]

As always, my pleasure! Thank you!

You can always ask a question! (With luck, I will always have an answer.)

The ‘Y_M’ calculation takes the known independent (‘data{1, i}(:,2)’) and dependent (‘data{1, i}(:,4)’) vector values and uses the newly-derived value for ‘X_M(i)’ to calculate (interpolate) the value for ‘Y_M(i)’. The only other additions were to define the figure object and the hold call before the loop, and plot the ‘(X_M(i),Y_M(i))’ values within the loop. The later xlim call simply ‘zooms’ the x-axis slightly to make the plotted ‘+’ markers more visible.