Ik was hier nog eens over na aan het denken en heb het eens gesimuleerd met behulp van Matlab met behulp van de formules van ZvdP van
bericht 54 van dit topic. Hierin wordt echter een verkeerde beginvoorwaarde gegeven:
\(\theta'(0) = \frac{-2 \pi}{24 \cdot 3600}\)
onafhankelijk van de aanvankelijke snelheid loodrecht op het aardoppervlak. Ik weet niet wat voor tijdstappen ik moet zetten, dus hier heb ik mee gespeeld. Zodra ik de tijdstappen te klein neem, wordt het systeem instabieler. Hij hangt dan bijvoorbeeld af van de beginhoek
\(\theta(0)\)
, dit kan natuurlijk niet.
Het resultaat wat ik vind, is dat wanneer je aanvankelijk vanaf de evenaar met 10 m/s omhoog springt, je een oostelijke afwijking krijgt van 1 mm en als je vanaf de evenaar vanaf 100 meter valt, je een westelijke afwijking krijgt van 22 mm.
Mijn vragen aan de lezer zijn:
* Klopt de code of heb ik ergens een fout gemaakt?
* Waarom wordt het al zo snel instabiel en hoe kan ik dit verbeteren?
De code die ik gebruikt heb staat hieronder afgebeeld. De berekeningen duren 1 seconde tot tijdstappen tot 1E-5 en 1 minuut tot tijdstappen tot 1E-7.
[attachment=6699:Jump_wit...imesteps.png]
[attachment=6700:Jump_wit...imesteps.png]
[attachment=6701:Drop_fro...imesteps.png]
[attachment=6702:Drop_fro...imesteps.png]
Code: Selecteer alles
% Set parameters
clc
clear all
R = 6.378E6; % Radius of earth in m
M = 5.9742E24; % Mass of earth in kg
G = 6.67300E-11; % Gravitational constant in m3 kg-1 s-2
r0 = R; % Initial r coordinate in m
v0 = 10; % Initial velocity upward in m s-1
expsteps = -2:-0.25:-4.5;
tsteps = 10.^expsteps;
T = zeros(length(tsteps),3);
theta0list = 0:0.5:3;
S = zeros(length(theta0list),length(tsteps));
for itheta0 = 1:length(theta0list);
theta0 = theta0list(itheta0);
for it = 1:length(tsteps)
tstep = tsteps(it);
% r = r coordinate of object
% rp = dr/dt (= r prime)
% rpp = d^2r/dt^2 (= r double prime)
% theta = theta coordinate of object
% thetap = dtheta/dt (= theta prime)
% thetapp = d^2theta/dt^2 (= theta double prime)
% Initial conditions
r = r0;
theta = theta0;
rp = v0;
thetap = -2*pi/(24*3600); % one round in 24 hours
tic; t = 0; % Reset the time
while r>=R
t = t+tstep; % Update the time
% Equation 1: m*M*G/r^2 = m(rpp - r*thetap^2)
rpp = -M*G/r^2 + r*thetap^2;
% Equation 2: 0 = 2*rp*thetap + r*thetapp
thetapp = -2*rp*thetap/r;
% Calculate first derivatives from second derivatives:
rp = rp + rpp*tstep;
thetap = thetap + thetapp*tstep;
% Calculate r and theta from first derivatives:
r = r + rp*tstep;
theta = theta + thetap*tstep;
end
thetaearth = theta0 + -2*pi/(24*3600)*t;
T(it,:) = [tstep (theta-thetaearth)*R toc];
end
disp(horzcat('Done in ',num2str(sum(T(:,3))),' seconds'))
S(itheta0,:) = T(:,2);
end
%% Plot the outcome
% Make the figure
close all
clear pm figprops
pm.ScreenSize = get(0,'ScreenSize');
pm.SizeRatio = 0.8;
figprops.Color = 'white';
figprops.Position = [pm.ScreenSize(3)*pm.SizeRatio * (1-pm.SizeRatio)/2 ...
pm.ScreenSize(3)*pm.SizeRatio * (1-pm.SizeRatio)/2 ...
pm.ScreenSize(3)*pm.SizeRatio ...
pm.ScreenSize(4)*pm.SizeRatio];
handles.figure = figure(figprops);
% Make the axes
clear axesprops
axesprops.LineWidth = 2;
axesprops.FontSize = 24;
axesprops.Parent = handles.figure;
axesprops.NextPlot = 'add';
axesprops.XScale = 'log';
axesprops.XDir = 'reverse';
handles.axes = axes(axesprops);
% Make the axis labels and title
handles.xaxis = get(handles.axes,'XLabel');
set(handles.xaxis,'FontSize',24);
set(handles.xaxis,'String','Time step (s)');
handles.yaxis = get(handles.axes,'YLabel');
set(handles.yaxis,'FontSize',24);
set(handles.yaxis,'String','Deviation to East (mm)');
handles.title = get(handles.axes,'Title');
set(handles.title,'FontSize',30);
set(handles.title,'String','Jump with 10 m/s');
% Make the plot
clear plotprops
plotprops.LineWidth = 2;
plotprops.Parent = handles.axes;
handles.plot = plot(T(:,1),1E3*S',plotprops);
for iplot = 1:length(handles.plot)
set(handles.plot(iplot),'DisplayName',...
horzcat('\theta_{0} = ',num2str(theta0list(iplot)), ' rad'));
end
legend(handles.plot)
pause; close all
Wat ik hoop is dat iemand anders de draad oppakt, en mijn afleidingen nakijkt. Als de afleidingen in de haak zijn, moet nog worden nagetrokken voor welke waarden van h en φB is voldaan aan de beperkende veronderstellingen die in de loop van de afleidingen onder meer ten aanzien van de grootte van de hoeken gemaakt zijn. De beperkende veronderstellingen staan in de afleidingen expliciet vermeld, en kom je bij het napluizen van de afleidingen daarom vanzelf tegen. Verder is het nog van belang dat wordt nagegaan hoe fouten in h en en φB in de uitkomst d doorwerken, en hoe groot het effect van de nauwkeurigheid van het gebruikte rekenprogramma op de nauwkeurigheid van de resultaten is. Ook zou men grafiekjes van d als functie van h voor verschillende breedtegraden φB kunnen maken. En men zou kunnen proberen deze functies in machtreeksen te ontwikkelen. Er is dus nog genoeg te doen! Ik lees graag in dit topic hoe dat uitpakt...
