Matlab code for method of 48 ordinates

Hi,

I have read about the

Method of 12 ordinates and the method of 24 ordinates

(Runge and Konig, Scarborough 1958)

Is there any matlab code that would work out 48 or 96 ordinates?

I have tried myself but it is very rough.

Thanks

ADT

Answer 1

[R]

Hi @Andrew Thorburn,

Here’s how you can implement the Method of 48 Ordinates in MATLAB using the function sin(x) and integrating it over the interval [0, 2pi].

function integral_value = ordinates_method(func, a, b, method_type)
    % func: the function to integrate
    % a: lower limit of integration
    % b: upper limit of integration
    % method_type: '48' or '96'

    % Define the number of ordinates based on the method type
    if strcmp(method_type, '48')
        n = 48;
        [ordinates, weights] = get_48_ordinates();
    elseif strcmp(method_type, '96')
        n = 96;
        [ordinates, weights] = get_96_ordinates();
    else
        error('Invalid method type. Use ''48'' or ''96''.');
    end

    % Map the ordinates to the interval [a, b]
    x = a + 0.5 * (b - a) * (ordinates + 1);  % Scale to [a, b]

    % Evaluate the function at the ordinates
    func_values = func(x);

    % Compute the integral using weighted sum
    integral_value = 0.5 * (b - a) * sum(weights .* func_values);
end

function [ordinates, weights] = get_48_ordinates()
    % Define 48 ordinates and their corresponding weights
    ordinates = linspace(-1, 1, 48);  % Example ordinates
    weights = ones(1, 48) * (2 / 48);  % Equal weights
end

function [ordinates, weights] = get_96_ordinates()
    % Define 96 ordinates and their corresponding weights
    ordinates = linspace(-1, 1, 96);  % Example ordinates
    weights = ones(1, 96) * (2 / 96);  % Equal weights
end

% Example usage: Integrate sin(x) from 0 to 2*pi using 48 ordinates
result = ordinates_method(@(x) sin(x), 0, 2*pi, '48');
disp(['Integral result of sin(x) from 0 to 2*pi: ', num2str(result)]);

Important considerations:

1. Function Definition: The ordinates_method function takes a function to integrate, the limits of integration, and the method type (either ‘48’ or ‘96’).
2. Ordinates and Weights: The get_48_ordinates and get_96_ordinates functions define the ordinates and weights. In this example, I’ve used simple linear spacing for the ordinates and equal weights for demonstration purposes.
3. Integration: The integral is computed by evaluating the function at the scaled ordinates and summing the weighted values.

[Andrew Thorburn]

Thanks very much, I'll look at that. I am interested in seeing how the coefficients change. Say a system

produces the same waveform time and time again. Then a sub system degrades, how do the coefficients change as the waveform slowly changes?

Answer 2

[Andrew Thorburn]

Hi, here is the code I wrote a few years ago now. I did try using the original paper written by Runge and Konig. I do not know German, so it was hard for me. I did get a result, but I had to a shift at the end. If someone could straighten it out for me that would be helpful. I am looking to pick it up again and run some experiments with changing values of components for a basic circuit. My coding skills are not the best:

%% *Ordinates calculations as per Runge and Konig*
% The x values are calculated by dividing 360 degrees by the size of the y data.
% Enter y values.
clear; clc;
y=[149 137 128 126 128 135 159 178 180 193 185 187 178 170 177 183 181 179 179 185 182 176 166 160 155 150 145 148 152 154 155 153 160 162 158 150 157 160 165 170 172 174 176 180 185 193 200 210];
sz = size(y);
r=sz(1,2);
x = [0:(360/r):360-(360/r)];
plot(x,y)

% Calculate a value for p. Number of data points divided by 4.
p = r/4;
% Calculate s and d.
s(1,1) = y(1,4*p);
s(1,2*p+1) = y(1,2*p);
for alphao = 1:1:2*p-1
    s(1,alphao+1) = y(1,alphao) + y(1,4*p-alphao);
    d(1,alphao) = y(1,alphao) - y(1,4*p-alphao);
end
%% 
% Moving onto the next calculations as defined in the next two images:
% 
% 
% 
% 
% 
% In the code I have changed the german letters to sh, dh, si, dl as shown 
% in the above image.
% 
% 

sh(1,p+1) = s(1,p+1);
for alphao = 1:1:p
    sh(1,alphao) = s(1,alphao) + s(1,(2*p+2)-alphao);
    dh(1,alphao) = s(1,alphao) - s(1,(2*p+2)-alphao);
end

si(1,p) = d(1,p);
for alphao = 1:p-1
    si(1,alphao) = d(1,alphao) + d(1,2*p-alphao);
    dl(1,alphao) = d(1,alphao) - d(1,2*p-alphao);
end
%% 
% Caculations can  now be done for the a cosine and b sine values:
% 
%  
% 
% Caculate the first line of (30c) when beta =0 and 2p. p=6 for these calculations.

suma=zeros(1,2*p+1);
sumb = zeros(1,2*p-1);

for alphao = 1:1:p+1
    suma(1,1) = suma(1,1) + sh(1,alphao)*cosd(0*x(1,alphao));
end
for alphao = 1:1:p+1
    suma(1,2*p+1) = suma(1,2*p+1) + sh(1,alphao)*cosd((2*p)*x(1,alphao));
end
%% 
% Calculate the a coefficients from second and fourth line of (30c) when 
% beta is even and odd = 1,2...2*p-1). Alpha runs from 0 to p, which is 1 to p+1 
% for the vector.

for beta = 2:1:2*p
    if rem(beta-1,2) == 0
        for alphao = 1:1:p+1 
            suma(1,beta) = suma(1,beta) + sh(1,alphao)*cosd((beta-1)*x(1,alphao));
        end
    else
        for alphao = 1:1:p 
            suma(1,beta) = suma(1,beta) + dh(1,alphao)*cosd((beta-1)*x(1,alphao));
        end
    end
end
%% 
% Calculate the b coefficients from the third and fifth lines of (30c) when 
% beta is even and odd and = 1,2,..2*p-1. 

for beta = 1:1:2*p-1
    if rem(beta,2) == 0
        for alphao = 1:1:p-1
            sumb(1,beta) = sumb(1,beta) + dl(1,alphao)*sind(beta*x(1,alphao+1));
        end
    else
        for alphao = 1:1:p 
            sumb(1,beta) = sumb(1,beta) + si(1,alphao)*sind(beta*x(1,alphao+1));
        end
    end
end

a(1,1) = suma(1,1)/r;
a(1,2*p+1) = suma(1,2*p+1)/r;
for alphao = 2:1:2*p
    a(1,alphao) = suma(1,alphao)/(r/2);
end
for alphao = 1:1:2*p-1
    b(1,alphao) = sumb(1,alphao)/(r/2);
end

ya = a(1,1)*ones(1, r);
for alphao = 2:1:2*p+1
    ya = ya + a(1,alphao)*cosd((alphao-1)*x);
end
yb = zeros(1, r);
for alphao = 1:1:2*p-1
    yb = yb + b(1,alphao)*sind(alphao*x);
end
yc=ya+yb;
yshift=circshift(yc,-1);
plot(x,yshift)
hold on
plot(x,y,'r')