I am computing the probablility of collision with the Monte Carlo method in MATLAB. The following is my code for it. The results don't seem to be accurate. Is there something wrong with my code? Also, it takes more than 5 min to run this code. Can you tell me why? %Initial Mean and Covariance mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000; P_01 = [6494.080 -376.139 0 0.0160 -0.494 0; -376.139 22.560 0 -9.883e-4 0.0286 0; 0 0 1.205 0 0 -6.071e-5; 0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0; -0.494 0.0286 0 -1.212e-6 3.762e-5 0; 0 0 -6.071e-5 0 0 3.390e-9]*1e-6; mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000; P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8; -376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9; -4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5; 0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13; -0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12; -5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]*1e-6; %Collision parameters R = 6378.0; %km mu = 398600.4415; %km^3/s^2 a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period1 = 2*pi*sqrt((a1^3)/mu); t_total1 = 0:10:2*period1; a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period2 = 2*pi*sqrt((a2^3)/mu); t_total2 = 0:10:2*period2; rho_AB = 15e-3; function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB) %Creating Eigen Matrix and Eigen Vector [R1,Lam1] = eig(P_01); [R2,Lam2] = eig(P_02); %Creating state vector x_01 = inv(R1)*(R1*mu_01 + sqrt(abs(diag(Lam1)))*randn(1,1000)); x_02 = inv(R2)*(R2*mu_02 + sqrt(abs(diag(Lam2)))*randn(1,1000)); x_final1 = zeros(6,1000); x_final2 = zeros(6,1000); % Propogate Orbit for i=1:1000 % Integrating the state using ode45 integrator init_x = x_01(:,i); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm1, init_x, options); x_final1(:,i) = x_out(length(tm1),:); end for j=1:1000 % Integrating the state using ode45 integrator init_x = x_02(:,j); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm2, init_x, options); x_final2(:,j) = x_out(length(tm2),:); end n_sampleA = 1000; n_sampleB = 1000; n_collide = 0; for idx1=1:n_sampleA for idx2=1:n_sampleB x1 = x_final1(1:3,idx1); x2 = x_final2(1:3,idx2); dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2); if dist <= rho_AB n_collide = n_collide + 1; end end end PC = n_collide/(n_sampleA*n_sampleB); % Plot fig1 = figure(); set(fig1, 'color', 'white'); scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:)); fig2 = figure(); set(fig2, 'color', 'white'); scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:)); end

I am computing the probablility of collision with the Monte Carlo method in MATLAB. The following is my code for it. The results don't seem to be accurate. Is there something wrong with my code? Also, it takes more than 5 min to run this code. Can you tell me why? %Initial Mean and Covariance mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000; P_01 = [6494.080 -376.139 0 0.0160 -0.494 0; -376.139 22.560 0 -9.883e-4 0.0286 0; 0 0 1.205 0 0 -6.071e-5; 0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0; -0.494 0.0286 0 -1.212e-6 3.762e-5 0; 0 0 -6.071e-5 0 0 3.390e-9]1e-6; mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000; P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8; -376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9; -4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5; 0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13; -0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12; -5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]1e-6; %Collision parameters R = 6378.0; %km mu = 398600.4415; %km^3/s^2 a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period1 = 2pisqrt((a1^3)/mu); t_total1 = 0:10:2period1; a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu)); period2 = 2pisqrt((a2^3)/mu); t_total2 = 0:10:2period2; rho_AB = 15e-3; function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB) %Creating Eigen Matrix and Eigen Vector [R1,Lam1] = eig(P_01); [R2,Lam2] = eig(P_02); %Creating state vector x_01 = inv(R1)(R1mu_01 + sqrt(abs(diag(Lam1)))randn(1,1000)); x_02 = inv(R2)(R2mu_02 + sqrt(abs(diag(Lam2)))randn(1,1000)); x_final1 = zeros(6,1000); x_final2 = zeros(6,1000); % Propogate Orbit for i=1:1000 % Integrating the state using ode45 integrator init_x = x_01(:,i); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm1, init_x, options); x_final1(:,i) = x_out(length(tm1),:); end for j=1:1000 % Integrating the state using ode45 integrator init_x = x_02(:,j); options = odeset('RelTol',1e-12, 'AbsTol',1e-12); [~, x_out] = ode45(@TwoBody, tm2, init_x, options); x_final2(:,j) = x_out(length(tm2),:); end n_sampleA = 1000; n_sampleB = 1000; n_collide = 0; for idx1=1:n_sampleA for idx2=1:n_sampleB x1 = x_final1(1:3,idx1); x2 = x_final2(1:3,idx2); dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2); if dist <= rho_AB n_collide = n_collide + 1; end end end PC = n_collide/(n_sampleA*n_sampleB); % Plot fig1 = figure(); set(fig1, 'color', 'white'); scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:)); fig2 = figure(); set(fig2, 'color', 'white'); scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:)); end

Principles of Heat Transfer (Activate Learning with these NEW titles from Engineering!)

8th Edition

ISBN:9781305387102

Author:Kreith, Frank; Manglik, Raj M.

Publisher:Kreith, Frank; Manglik, Raj M.

Chapter5: Analysis Of Convection Heat Transfer

Section: Chapter Questions

Problem 5.19P: 5.19 Suppose that the graph below shows measured values of for air in forced convection over a...

See similar textbooks

Related questions

Question

%Initial Mean and Covariance

mu_01 = [153446.180; 41874155.872;0;3066.875;-11.374;0]/1000;

P_01 = [6494.080 -376.139 0 0.0160 -0.494 0;

-376.139 22.560 0 -9.883e-4 0.0286 0;

0 0 1.205 0 0 -6.071e-5;

0.0160 -9.883e-4 0 4.437e-8 -1.212e-6 0;

-0.494 0.0286 0 -1.212e-6 3.762e-5 0;

0 0 -6.071e-5 0 0 3.390e-9]*1e-6;

mu_02 = [153446.679; 41874156.372;5.000;3066.865;-11.374;-1.358e-6]/1000;

P_02 = [6494.224 -376.156 -4.492e-5 0.0160 -0.494 -5.902e-8;

-376.156 22.561 2.550e-6 -9.885e-3 0.0286 3.419e-9;

-4.492e-5 2.550e-6 1.205 -1.180e-10 3.419e-9 -6.072e-5;

0.0160 -9.883e-3 -1.180e-10 4.438e-8 -1.212e-6 -1.448e-13;

-0.494 0.0286 3.419e-9 -1.212e-6 3.762e-5 4.492e-12;

-5.902e-8 3.419e-9 -6.072e-5 -1.448e-13 4.492e-12 3.392e-9]*1e-6;

%Collision parameters

R = 6378.0; %km

mu = 398600.4415; %km^3/s^2

a1 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu));

period1 = 2*pi*sqrt((a1^3)/mu);

t_total1 = 0:10:2*period1;

a2 = 1/((2/norm(mu_01(1:3)))-((norm(mu_01(4:6))^2)/mu));

period2 = 2*pi*sqrt((a2^3)/mu);

t_total2 = 0:10:2*period2;

rho_AB = 15e-3;

function [PC, n_collide] = MonteCarlo_PC(tm1,P_01,mu_01,tm2,P_02,mu_02,rho_AB)

%Creating Eigen Matrix and Eigen Vector

[R1,Lam1] = eig(P_01);

[R2,Lam2] = eig(P_02);

%Creating state vector

x_01 = inv(R1)*(R1*mu_01 + sqrt(abs(diag(Lam1)))*randn(1,1000));

x_02 = inv(R2)*(R2*mu_02 + sqrt(abs(diag(Lam2)))*randn(1,1000));

x_final1 = zeros(6,1000);

x_final2 = zeros(6,1000);

% Propogate Orbit

for i=1:1000

% Integrating the state using ode45 integrator

init_x = x_01(:,i);

options = odeset('RelTol',1e-12, 'AbsTol',1e-12);

[~, x_out] = ode45(@TwoBody, tm1, init_x, options);

x_final1(:,i) = x_out(length(tm1),:);

end

for j=1:1000

% Integrating the state using ode45 integrator

init_x = x_02(:,j);

options = odeset('RelTol',1e-12, 'AbsTol',1e-12);

[~, x_out] = ode45(@TwoBody, tm2, init_x, options);

x_final2(:,j) = x_out(length(tm2),:);

end

n_sampleA = 1000;

n_sampleB = 1000;

n_collide = 0;

for idx1=1:n_sampleA

for idx2=1:n_sampleB

x1 = x_final1(1:3,idx1);

x2 = x_final2(1:3,idx2);

dist = sqrt((x1(1)-x2(1))^2 + (x1(2)-x2(2))^2 + (x1(3)-x2(3))^2);

if dist <= rho_AB

n_collide = n_collide + 1;

end

PC = n_collide/(n_sampleA*n_sampleB);

% Plot

fig1 = figure();

set(fig1, 'color', 'white');

scatter3(x_final1(1,:),x_final1(2,:),x_final1(3,:));

fig2 = figure();

set(fig2, 'color', 'white');

scatter3(x_final2(1,:),x_final2(2,:),x_final2(3,:));

end

Expert Solution

Step by step

Solved in 2 steps

SEE SOLUTION Check out a sample Q&A here