Including a periodic piecewise function of time in coupled ODE
Ältere Kommentare anzeigen
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.
and are over the same domain/variable x.Thank you for any help,
Tom
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
Thomas Dixon
am 2 Feb. 2021
Thomas Dixon
am 2 Feb. 2021
Alan Stevens
am 2 Feb. 2021
What you want should work ok if you 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)];
[x,A] = ode45(dA, [0 3*m] ,[Ap0; As0; 0; 0]);
However, ode45 can't meet the integration tolerances at around x = 353. This suggests you should take a very careful look at your (complicated) function, and/or your data values.
Thomas Dixon
am 2 Feb. 2021
Thomas Dixon
am 2 Feb. 2021
Walter Roberson
am 5 Feb. 2021
piecewise functions can only be used with ode*() functions if they have C2 or higher continuity (their second derivatives must be continuous.) The ode*() routines will typically not detect the problem if the functions are C1 continuous (first derivative is continuous) and will instead silently produce the wrong value. The ode*() routines will typically not detect problems with C1 continuity if the jump is proportionately less than the current integration tolerance, and again will silently return the wrong value in such a case.
Walter Roberson
am 5 Feb. 2021
Ah the x vector is spaced non-linearly. I have not seen this before. Is this something unphysical?
The ode*() functions do not return values at constant time (or in your case, x) intervals unless you pass in a tspan that has more than 2 elements and the tspan is spaced regularly.
The ode*() functions are all adaptive functions, changing the step interval as needed to meet local conditions. In places where the integration is going well, it will increase the step size. It will keep increasing the step size as long as it is successful in doing so. The functions make calculations up to a certain order, and use the information to deduce a new state. They then use that deduced state to predict at one or more additional locations, and evaluate at those locations and cross-check the actual values against predictions; if the result is within tolerance then the deduced state is "accepted" (and next round, the step size will be increased.) If the actual values were not within tolerance of those predicted by the deduced state, then the trial step is rejected, and the step size is decreased and a new series of measurements is taken relative to the previous state. So as long as everything goes well, the step size accelerates, and when it figures out that it has overstepped, it retreats back to a smaller step.
This is why the x vector is spaced non-linearly: the spacing used is adaptive, more sparse when that works, less sparse when the conditions are changing rapidly.
Thomas Dixon
am 5 Feb. 2021
Walter Roberson
am 5 Feb. 2021
In particular you would need to match the derivatives at the boundaries.
If the derivatives at the boundaries are matched, in some cases you can use argument reduction, if the periods are all the same length and are identical. If, though, the duration stretches as you go, or the amplitude changes as you go, argument reduction alone might not work.
Linear interpolation is not C1 continuous; you need spline for C1 and C2 continuity.
Thomas Dixon
am 5 Feb. 2021
Bearbeitet: Thomas Dixon
am 5 Feb. 2021
Walter Roberson
am 5 Feb. 2021
boundaries of each piecewise section.
Thomas Dixon
am 5 Feb. 2021
I am very grateful for the depth of this, I think I may have confused the problem.
I want to solve some coupled differential equation of the form:
dA(1)/dx = -eta[x]*A(2)[x]
dA(2)/dx = eta[x]*A(1)[x]
The integration function I want to find is A, and while eta[x] is a function of x I wouldn't care if the value given to the ODE solver is just a value of eta defined at certain x points (these points are then defined by where the ODE solver is numerically computing the derivative), then the ODE solver treats eta as just a number which it may/may not need to find depending on where it wants to solve the ODE. (I hope that makes sense, I do understand that a constant that changes is not a constant, but I am looking for insight into how this may be done without defining it in some crude made up way from my head)
But to be clear the the function eta is meant to be continuous, differentiable and single valued, as shown in the plot:
I believe it is very similar to a two sided bump function which has no analytic form for its entire period (I am happy to be proven wrong). I may then want to extend this (see eta full) which would then definitely not have an analytic form (Very happy to be proven wrong).
The ODE solver works very well in the period of eta(x) this is because eta(x) is defined in that period. However then eta(x) goes to 0 because it hasnt been repeated (that is eta_full). but now eta_full is not a function of x and so does not go into the ODE solver. The reason the ODE solver can comlete solving though is because MATLAB has defined the peicewise function eta to go to zero outside of its period. This makes solving the differential equation very easy as now al derivatives are proportional to zero, the function doesnt change and you can see the reuslt in this plot:
Where the blue circles are solutions when eta is not included and the coloured lines are when eta is included. Note the integrator did not 'fail' it was the combination of an ill defined function as a coefficient, an illdefined ODE and my superior lack of pragramming ability, I believe this can be sorted but I am not making my intentions clear enough.
I appreciate your continued support,
Tom
Ps this is obviously going to be extended so I would like a solution where I don't have to define the peicewise function for the entire duration by stitching it together multuiple times so some repeated functionality would be desirable.
Walter Roberson
am 5 Feb. 2021
Argument reduction.
eta(ReduceEtaArgument(x))
where ReduceEtaArgument(x) is something like
mod(x, duration_of_one_eta_cycle)
Alan Stevens
am 5 Feb. 2021
You can repeat your eta function shape for as long as you like using the following function
function etaval = eta(x)
a=0.15;
b=0.2;
m = 448;
x = mod(x,m);
etaval= (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))));
end
Your dA function still gives problems around x = 352 though.
Thomas Dixon
am 5 Feb. 2021
Bearbeitet: Thomas Dixon
am 5 Feb. 2021
Thomas Dixon
am 5 Feb. 2021
Aplogies here are the plots
Thomas Dixon
am 5 Feb. 2021
Weitere Antworten (0)
Kategorien
Mehr zu Loops and Conditional Statements finden Sie in Hilfe-Center und File Exchange
am 1 Feb. 2021
am 5 Feb. 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
