Funktionen mit Matlab < Matlab < Mathe-Software < Mathe < Vorhilfe
|
| Status: |
(Frage) überfällig  | | Datum: | 16:26 Di 10.03.2009 | | Autor: | Mirco |
Hallo, ich möchte gerne einen Matlabcode schreiben, komme allerdings mathematisch nicht weiter. Ich möchte eine Rechteckwelle simulieren, die aus 2 Schwingungen besteht, die jedoch verschieden hohe x-Werte haben. Außerdem weiss ich nicht, wie ich diese Schwingungen ungleichmäßig mache - sprich wie ich eine große Rechtwelle (weiter Bereich auf der x-Achse) und eine kleine Rechteckwelle (kleiner Bereich auf der x-Achse) mathematisch definiere. Bislang habe ich die Rechteckwelle nur mit gleichgroßen x-Werten geschafft.
Der MatLab-Code für die bisherige Variante sieht so aus:
%%% To run: erode( flag, amplitude, frequency, pedestal )
%%% where, flag defines erosion rate:
%%% 1 = zero (steady state)
%%% 2 = constant
%%% 3 = full-wave rectification
%%% 4 = square-wave
%%% 5 = saw-tooth
%%% 6 = exponential
%%% 7 = step function
%%% and amplitude/frequency/pedastal control the periodic erorison
%%% sequences (3-7) where erosion units are in m/MYrs. Amplitude
%%% defines max erosion and pedastal defines min erosion
%%%
function[concentration] = erode( flag, AMPLITUDE, FREQUENCY, PEDESTAL, ATLANTICBOREAL, ANTHROPOGEN, RECENT );
%%% DEFAULT PARAMETERS OF PERIODIC EROSIONS
%%%
if( ~exist( 'AMPLITUDE' ) ) AMPLITUDE = 150; end;
if( ~exist( 'FREQUENCY' ) ) FREQUENCY = 1; end;
if( ~exist( 'PEDESTAL' ) ) PEDESTAL = 22; end;
if( ~exist( 'ATLANTICBOREAL' ) ) ATLANTICBOREAL = 22; end;
if( ~exist( 'ANTHROPOGEN' ) ) ANTHROPOGEN = 150; end;
if( ~exist( 'RECENT' ) ) RECENT = 22; end;
%%% VARIOUS CONSTANTS
%%%
flag = 4;
timestep = 10; % time step - this times the below defines the timestep
% for high erosion rates, lower the time step
maxtime = 10e4; % max time - sets the run duration of the model
time = [timestep : timestep : maxtime]; % time interval
P0 = 17*timestep; % surface production
C0 = 100; % initial concentration
mu = 2.65/165; % constant
lambda = log(2)/1.5e6 * timestep;% decay rate
e = erosionrate( flag, time, AMPLITUDE, FREQUENCY, PEDESTAL, ATLANTICBOREAL, ANTHROPOGEN, RECENT );
e = 1e2/1e6 * e; % m/Myrs to cm/yrs
e = e * timestep;
%%% ERODE IT...
%%%
%%%100 sets the sample interval for the data - reduce for finer resolution
ttt = 1 : 1 : length(time);
concentration = zeros( 1, length(ttt) );
count = 1;
for t = ttt
C = C0;
CS = fliplr( cumsum( fliplr(e(1:t)) ) );
for k = 1 : t
c = P0 * exp(-mu*CS(k)) - lambda*C;
C = C + c;
end
concentration(count) = C;
count = count + 1;
end
%%% ESTIMATE EROSION FROM ASSUMED STEADY-STATE CONCENTRATION
%%%
eEst = P0 ./ (concentration*mu) -lambda;
eEst = eEst/timestep * 1e6/1e2;
%%% PLOT EVERYTHING
%%%
set( gca, 'FontSize', 12 );
subplot(211); plot( time(ttt), concentration, 'bo-' );
xlabel( 'Time (y)' );
ylabel( '^{10}Be (atoms/gram)' );
set( gca, 'xlim', [1e3 maxtime]);
set( gca, 'ylim', [10 2e5] ); %change this to reflect expected max concentration
grid on;
subplot(212); plot( time(ttt), e(ttt)/timestep*1e6/1e2, 'r' );
hold on;
plot( time(ttt), eEst, 'k--' );
hold off;
%set( gca, 'Xscale', 'log' );
xlabel( 'Time (y)' );
ylabel( 'Erosion (m/My)' );
set( gca, 'xlim', [1e3 maxtime] ); %change this to 1e5 when time is [mm] >10^6 [/mm] yrs-1e2 otherwise
set( gca, 'ylim', [10 4e2] ); %change this to be just over the max erosion
format long;
disp( [time(ttt)' concentration' e(ttt)'/timestep*1e6/1e2 eEst'] );
return;
%%% =======================================================================
%%% DEFINE ERORSION RATE AS A FUNCTION OF TIME
%%%
function[e] = erosionrate( flag, time, amp, freq, ped, atl, ant, rec )
if( flag == 1 ); % zero (steady state)
e = zeros(size(time));
elseif( flag == 2 ) % constant
e = amp * ones(size(time));
elseif( flag == 3 ) % full-wave rectification
e = sin( freq*2*pi*time/max(time) );
e(find(e<0)) = 0;
e = (amp-ped) * e + ped;
elseif( flag == 4 ) % square-wave
e = sin( freq/1.79*2*pi*time/max(time) );
e( find(e<0) ) = 0;
e( find(e>0) ) = 1;
e = (amp-ped) * e + ped;
elseif( flag == 5 ) % saw-tooth
N = length(time);
n = N / freq;
saw = [1:n]/n;
e = repmat( saw, 1, freq );
e = (amp-ped) * e + ped;
if( length(e)<N ) % freq is a non-integer multiple of time
e( length(e)+1:N ) = e(1:N-length(e));
end
elseif( flag == 6 ) % exponential
N = length(time);
n = N / freq;
expo = exp(-([n:-1:1]/n)*6);
e = repmat( expo, 1, freq );
e = (amp-ped) * e + ped;
if( length(e)<N ) % freq is a non-integer multiple of time
e( length(e)+1:N ) = e(1:N-length(e));
end
elseif( flag == 7 ) % step
N = length(time);
n = round(N/1.5);
e = zeros( 1,n );
e( length(e)+1:N ) = 1;
e = (ped-amp) * e + amp;
elseif( flag == 8 ) % constant
e = amp * ones(size(time));
else
fprintf( 'erosionrate(): invalid flag (%d)-using steady [mm] state\n', [/mm] flag );
e = zeros(size(time));
end
flag bezieht sich auf eine alte Auswhl, ist nicht mehr wichtig!
Die alte Rechteckwelle (rot in der Darstellung) wurde über amp, ped, und freq definiert. Ich habe nun die neuen Definitionen atl, anth, und rec eingefügt die die weiterführende Rechteckwelle definieren sollen, was ich aber nicht hinbekommen habe.
Im Endeffekt soll die Rechteckwelle folgendes Aussehen haben: Werte auf der y-Achse (x-Achsen-Werte in Klammern): 150 (0 bis 8,95) - 22 (8,95 bis 9,65) - 120 ( 9,65 bis 9,9925) - 22 (9,9925 bis 10)!
Wenn mir jemand helfen kann bin ich Ihm auf Ewig dankbar!!!
Gruss Mirco
Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 17:20 Mi 25.03.2009 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|
