Including a periodic piecewise function of time in coupled ODE
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Thomas Dixon
am 1 Feb. 2021
Kommentiert: Thomas Dixon
am 5 Feb. 2021
Hi All,
I want to solve an ODE which includes a function which is periodic and peicewise contructed (although not with the peicewise functionality...). This is something very similar to:
However not the same, I think.
I have made a go at trying to have the function defined symbolically with in its period and as just a general function for a number of periods. I then show how I solve the coupled ODE's without this function as a coefficient.
I have then commented out the code which is my (wrong) way of trying to incorporate the function into the ODE45 so that the code will run.
At the end I have plotted my solution for the working ODE without the function and the function itself to demonstrate that both are solved over the same x-domain.
The peicewise function is in black, and as you can see it is smooth, continuoud and differntiable.
I am not adverse to redefining the function if it can be done, or smoothing or interpolating the function in the ODE solver.
clear
a=0.15;
b=0.2;
m=448;
x=linspace(0,m,10000);
intvl = [0 3*m];
eta= @(x) [(0<=x & x<m/2).*(-a.*log(exp(-(2.*x)./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-(2.*x))./(a.*b.*m)))) + (m/2<=x & x<m).*(a.*log(exp(-(2*(x-m/2))./(a.*b.*m))+exp(-(1)./(a))+exp(-(m-2.*(x-m/2))./(a.*b.*m))))];
etafull = repmat(eta(x),1,3);
xfull=linspace(intvl(1), intvl(2),length(etafull));
figure(1)
plot(xfull,etafull)
xlim([0 2.5*m])
w0 = 2*pi*67*10^9;
wj = 2*pi*86*10^9;
wp = 2*pi*12*10^9;
wsfac = 0.6;
wifac = 1-wsfac;
ws = wsfac.*wp;
w2p = 2*wp;
Ap0 = 0.5*w0/wp;
As0 = Ap0*sqrt(0.0057*wp/ws);
%As0
A2p0 = 0;
wi = wifac.*wp;
kp = (wp/w0)*(1/(sqrt(1-(wp/wj)^2)));
ks = (ws/w0)*(1/(sqrt(1-(ws/wj)^2)));
ki = (wi/w0)*(1/(sqrt(1-(wi/wj)^2)));
k2p = (w2p/w0)*(1/(sqrt(1-(w2p/wj)^2)));
delk = 3*ws*wi*wp/(2*w0*(wj^2));
modk = sqrt(wp.*Ap0^2/(ws*As0^2 + wp.*Ap0^2));
maxBeta=0.433
dA = @(xfull,A)[-(maxBeta/2)*ks*ki*A(2)*A(3)*exp(1i*(ks+ki-kp)*xfull) + (maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*xfull);
(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*(kp-ki-ks)*xfull);
(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*(kp-ks-ki)*xfull);
-(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*xfull)];
[xfull,A] = ode45(dA, xfull ,[Ap0; As0; 0; 0]);
% dA = @(x,A)[-eta*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
% eta*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
% eta*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
% -eta*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];
%
% [x,A] = ode45(dA, x ,[Ap0; As0; 0; 0]);
P=[(wp.^2)*(abs(A(:,1))).^2/((wp.^2)*(abs(A(1,1))).^2), (ws.^2)*(abs(A(:,2)).^2)/((wp.^2)*(abs(A(1,1))).^2), (wi.^2)*(abs(A(:,3)).^2)/((wp.^2)*A(1,1).^2), (w2p.^2)*(abs(A(:,4)).^2)/((wp.^2)*A(1,1).^2)];
figure(2)
plot(xfull,P)
hold on
plot(xfull,etafull,'k','LineWidth',3)
hold off
% dA = @(x,A,eta)[-eta*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
% eta*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
% eta*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
% -eta*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];
%
% [x,A] = ode45(dA, x ,[Ap0; As0; 0; 0]);
Post edit: the variable x is the same variable to solve for A and that eta is defined over, that is and are over the same domain/variable x.
Thank you for any help,
Tom
0 Kommentare
Akzeptierte Antwort
Alan Stevens
am 1 Feb. 2021
Because eta is itself a function of x, you need to have
dA = @(x,A)[-eta(x)*(maxBeta/2)*ks*ki*A(2)*A(3)*exp(-1i*delk*x) + ((k2p-kp)/(kp*(1-(wp/wj)^2)))*(maxBeta/2)*k2p*A(4)*kp*conj(A(1))*exp(1i*(k2p-2*kp)*x);
eta(x)*(maxBeta/2)*ki*kp*conj(A(3))*A(1)*exp(1i*delk*x);
eta(x)*(maxBeta/2)*ks*kp*conj(A(2))*A(1)*exp(1i*delk*x);
-eta(x)*(maxBeta/4)*kp^2*A(1)^2*exp(1i*(2*kp-k2p)*x)];
17 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!