Clear Filters
Clear Filters

Numerically integrating a symbolic equation in a for loop?

5 views (last 30 days)
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

Accepted Answer

Harvey Rael
Harvey Rael on 25 Jun 2018
Fixed it - simply needed to flick the 'ArrayValued' setting in integral to true.
  1 Comment
Torsten
Torsten on 26 Jun 2018
I don't see the reason why this should be necessary in your code.
Best wishes
Torsten.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 25 Jun 2018
In the definition of "R" and "thetadot", you use "k1" instead of "k1(i)".
  6 Comments
Harvey Rael
Harvey Rael on 25 Jun 2018
Post your code then? Is y a vector or a single point solution? Because I contiually get errors
Torsten
Torsten on 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.

Sign in to comment.

Categories

Find more on Mathematics in Help Center and File Exchange

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!