Runge-Kutta 4th order method for 2nd order differential equation with complex number
15 次查看(过去 30 天)
显示 更早的评论
Hi,
I want to solve following equation, which has complex part:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546607/image.png)
I tried Runge-Kutta 4th order method but it does not provide me with correct result.
I assume:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546612/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546617/image.png)
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):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1546627/image.jpeg)
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 个评论
回答(1 个)
Fabio Freschi
2023-11-22
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 个评论
Fabio Freschi
2023-11-23
If you want you can share the original document with the formulas... fresh eyes can help
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!