Including a periodic piecewise function of time in coupled ODE

2 ビュー (過去 30 日間)
Thomas Dixon
Thomas Dixon 2021 年 2 月 1 日
コメント済み: Thomas Dixon 2021 年 2 月 5 日
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.
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));
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);
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));
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);
[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)];
hold on
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,


Alan Stevens
Alan Stevens 2021 年 2 月 1 日
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);
  17 件のコメント
Thomas Dixon
Thomas Dixon 2021 年 2 月 5 日
Ignore that last bit it is a fundemental misunderstanding of what a function and input arguments are.
I guess the ODE solver also goes to this function at any x point it wants and asks 'what the value is' arbitrarily.
It is then simple I assume to pass a vectorised x (ie the one each ode solver returns me x1, x2) to the same function and simply plot the result that is:
eta(0,m/2,m) = [0,0,0]
eta(0:1:m/2) = ['whatever numbers make that tapered sine wave evaluated at x=0:1:m/2']
so I can always ask for
and then:
to see eta.


その他の回答 (0 件)


Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by