Runge-Kutta 4th order method for 2nd order differential equation with complex number
19 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi,
I want to solve following equation, which has complex part:
I tried Runge-Kutta 4th order method but it does not provide me with correct result.
I assume:
I have used following code:
close all; clear variables;
h=0.01e-6; % step size
x = (-5e-6):h:(20e-6);
y = zeros(1,length(x));
z = zeros(1,length(x));
y(1) = 10; % y at x=-5e-6
z(1) = -1e-1; % dy/dx at x=-5e-6
dWdx = @(x) 2.3139e+13 * exp((-x.^2)/(2e-12));
W = @(x) integral(dWdx,-100e-6,x);
A = @(x) 7.8957e+03 *W(x);
B = @(x) 1/W(x)*dWdx(x);
f = @(x,y,z) z;
g = @(x,y,z) B(x)*z+ 1i*A(x)*y;
N=length(x);
for j=1:(N-1)
k0 = h*f(x(j),y(j),z(j));
L0 = h*g(x(j),y(j),z(j));
k1= h*f(x(j)+0.5*h,y(j)+0.5*k0,z(j)+0.5*L0);
L1= h*g(x(j)+0.5*h,y(j)+0.5*k0,z(j)+0.5*L0);
k2= h*f(x(j)+0.5*h,y(j)+0.5*k1,z(j)+0.5*L1);
L2= h*g(x(j)+0.5*h,y(j)+0.5*k1,z(j)+0.5*L1);
k3= h*f(x(j)+1*h,y(j)+1*k2,z(j)+1*L2);
L3= h*g(x(j)+1*h,y(j)+1*k2,z(j)+1*L2);
y(j+1)=y(j)+1/6*(k0+2*k1+2*k2+k3);
z(j+1)=z(j)+1/6*(L0+2*L1+2*L2+L3);
end
figure();
% plot(x,y./max(y));
plot(x,y);
xlim([-5 20]*1e-6);
grid on; hold on;
The solution should have following shape (see normalized red line, for 1GHz):
Can someone help me with this scenario? There where code (based on which i had started) that solves 2nd order ODE with 4rd order Runge-Kutta (https://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode), however I couldnt find a solution for equation with complex number.
Thanks in advance :)
0 commentaires
Réponses (1)
Fabio Freschi
le 22 Nov 2023
Your code is working correctly: in fact, if you use the built-in ode45 solver, you obtain exactly the same result. My guess there is a mistake in the definition of your coefficients W, A, B, or the initial conditions. Note that if you are not asked to write your own code, ode45 is more efficient and accurate, since it is using adaptive step size with error control.
close all; clear variables;
h=0.01e-6; % step size
x = (-5e-6):h:(20e-6);
y = zeros(1,length(x));
z = zeros(1,length(x));
y(1) = 10; % y at x=-5e-6
z(1) = -1e-1; % dy/dx at x=-5e-6
dWdx = @(x) 2.3139e+13 * exp((-x.^2)/(2e-12));
W = @(x) integral(dWdx,-100e-6,x);
A = @(x) 7.8957e+03 *W(x);
B = @(x) 1/W(x)*dWdx(x);
% use built in ode45
dfdx = @(x,f)[B(x)*f(1)+ 1i*A(x)*f(2); f(1)];
[t,f] = ode45(dfdx,[x(1) x(end)],[z(1) y(1)]);
% plot
h0 = figure;
plot(t,f(:,2));
xlim([-5 20]*1e-6);
grid on; hold on;
% OP code to compare
f = @(x,y,z) z;
g = @(x,y,z) B(x)*z+ 1i*A(x)*y;
N=length(x);
for j=1:(N-1)
k0 = h*f(x(j),y(j),z(j));
L0 = h*g(x(j),y(j),z(j));
k1= h*f(x(j)+0.5*h,y(j)+0.5*k0,z(j)+0.5*L0);
L1= h*g(x(j)+0.5*h,y(j)+0.5*k0,z(j)+0.5*L0);
k2= h*f(x(j)+0.5*h,y(j)+0.5*k1,z(j)+0.5*L1);
L2= h*g(x(j)+0.5*h,y(j)+0.5*k1,z(j)+0.5*L1);
k3= h*f(x(j)+1*h,y(j)+1*k2,z(j)+1*L2);
L3= h*g(x(j)+1*h,y(j)+1*k2,z(j)+1*L2);
y(j+1)=y(j)+1/6*(k0+2*k1+2*k2+k3);
z(j+1)=z(j)+1/6*(L0+2*L1+2*L2+L3);
end
figure(h0), hold on;
% plot(x,y./max(y));
plot(x,y);
xlim([-5 20]*1e-6);
grid on; hold on;
6 commentaires
Fabio Freschi
le 23 Nov 2023
If you want you can share the original document with the formulas... fresh eyes can help
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!