MATLAB Answers

Including a periodic piecewise function of time in coupled ODE

62 views (last 30 days)
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

Accepted Answer

Alan Stevens
Alan Stevens on 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 Comments
Thomas Dixon
Thomas Dixon on 5 Feb 2021
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
eta(some_vector)
and then:
plot(some_vector,eta)
to see eta.
Lovely

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!