Numerically integrating a symbolic equation in a for loop?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Harvey Rael
am 25 Jun. 2018
Kommentiert: Torsten
am 26 Jun. 2018
Hi! I'll try to explain qualitatively first what I am trying to do before then giving code. I have an equation that depends on a variable (theta) that changes (though how it changes doesn't matter). However, theta's end point depends on another variable N, and N is significant. I have a function that needs to be numerically integrated in terms of theta, to its end point which is pi/2+2piN. I want to change N with a for loop and store the outputs to be plotted.
Here's my code:
clc;clear;
% Define constants
AU = astroConstants(2);
mu = astroConstants(4);
% Define start and finish parameters for the exponential sinusoid.
r1 = AU; % Initial radius
psi = pi/2; % Final polar angle of Mars/finish transfer
phi = pi/2;
r2 = 1.5*AU; %Final radius
k2=0.2;
N = 0:1:5;
k1 = zeros(size(N));
y = zeros(size(N));
for i = 1:length(N)
k1(i) = sqrt( ( (log(r1/r2) + sin(k2*(psi + 2*pi.*N(i)))*tan(0)/k2) / (1-cos(k2*(psi+2*pi.*N(i)))) )^2 +
tan(0)^2/k2^2 );
k0 = r1./exp(k1(i)*sin(phi)); %The parameters of the exponential sinusoid
R = @(theta) k0*exp(k1(i)*sin(k2.*theta + phi)); %Equation of the exponential sinusoid
theta_dot = @(theta) sqrt((mu./(R(theta).^3))./((tan(0))^2 + k1(i)*k2^2*sin(k2*theta + phi) + 1));
z = @(theta)1./theta_dot(theta); %change of variable name for ease of use
y(i) = integral(z,0,(psi+2.*pi.*N(i))) - 94608000; %function to be plotted
end
The error message I keep getting is :
"Error in integral (line 88) Q = integralCalc(fun,a,b,opstruct);
Error in Untitled2 (line 23) y(i) = integral(theta_dot,0,(psi+2.*pi.*N(i))) - 94608000; %function to be plotted"
Anyone know where I am going wrong? When I have N as a single number, everything works fine.
EDIT : As @Torsten pointed out I forgot add in a couple of (i)'s in my code - so I have added them in. Error message persists
0 Kommentare
Akzeptierte Antwort
Harvey Rael
am 25 Jun. 2018
1 Kommentar
Torsten
am 26 Jun. 2018
I don't see the reason why this should be necessary in your code.
Best wishes
Torsten.
Weitere Antworten (1)
Torsten
am 25 Jun. 2018
In the definition of "R" and "thetadot", you use "k1" instead of "k1(i)".
6 Kommentare
Torsten
am 26 Jun. 2018
function main
AU = 2.0;
mu = 0.5;
% Define start and finish parameters for the exponential sinusoid.
r1 = AU; % Initial radius
psi = pi/2; % Final polar angle of Mars/finish transfer
phi = pi/2;
r2 = 1.5*AU; %Final radius
k2 = 0.2;
N = 0:1:5;
k1 = zeros(size(N));
y = zeros(size(N));
for i = 1:length(N)
k1(i) = sqrt( ( (log(r1/r2) + sin(k2*(psi + 2*pi.*N(i)))*tan(0)/k2) / (1-cos(k2*(psi+2*pi.*N(i)))) )^2 +...
tan(0)^2/k2^2 );
k0 = r1./exp(k1(i)*sin(phi)); %The parameters of the exponential sinusoid
R = @(theta) k0*exp(k1(i)*sin(k2.*theta + phi)); %Equation of the exponential sinusoid
theta_dot = @(theta) sqrt((mu./(R(theta).^3))./((tan(0))^2 + k1(i)*k2^2*sin(k2*theta + phi) + 1));
z = @(theta)1./theta_dot(theta); %change of variable name for ease of use
y(i) = integral(z,0,(psi+2.*pi.*N(i))) - 94608000; %function to be plotted
end
end
Best wishes
Torsten.
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!