Numerically integrating a symbolic equation in a for loop?

5 views (last 30 days)
Harvey Rael on 25 Jun 2018
Commented: Torsten on 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;
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

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

Torsten on 25 Jun 2018
In the definition of "R" and "thetadot", you use "k1" instead of "k1(i)".
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 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;
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.

Categories

Find more on Mathematics in Help Center and File Exchange

R2018a

Community Treasure Hunt

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

Start Hunting!