62 views (last 30 days)

Show older comments

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

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)];

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

Start Hunting!