Matlab code for method of 48 ordinates

43 次查看（过去 30 天）
Andrew Thorburn
Andrew Thorburn 2024-10-25，9:13
编辑： Voss 2024-10-25，18:14
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

请先登录，再回答此问题。

采纳的回答

R
R 2024-10-25，11:12
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.
  1 个评论
Andrew Thorburn
Andrew Thorburn 2024-10-25，12:28
编辑：Torsten 2024-10-25，17:59
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?

请先登录，再进行评论。

更多回答（1 个）

Andrew Thorburn
Andrew Thorburn 2024-10-25，14:22
编辑：Voss 2024-10-25，18:14
Ran in:
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')

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by