Solving the differential equation by the Runge-Kutta method (ode45)

4 次查看（过去 30 天）
Babr
Babr 2024-3-6
评论： Babr 2024-3-6
I want solve this equation numerically and use fourth order Runge-Kutta method. (ode45)
What is the code?
My code:
function dfdt = odefun2(x,y)
dfdt = [y(2); y(1) .^(-3) - (w .*r0 ./c) .^2 .*phi2 .*y(1)];
end
****
clc, clear, close all;
w = 2e16;
r0 = 30e-6;
c = 3e8;
p = 2.5;
f = 1;
Wp0 = p .*c ./r0;
a02 = 0.050;
phi2 = (Wp0 ./w) .^2 .*(3 .*a02 ./f .^4) ./4 ./(1 + a02 ./2 ./f .^2) .^1.5...
.*(1 + (56 + a02 ./f .^2) ./3 ./(1 + a02 ./2 ./f .^2) ./p .^2 ./f .^2);
[t, y] = ode45(@odefun2,[0 3],[1,0]);
******
But it gives an error ..
I have attached two pictures to the question, look at them if you need.
  5 个评论
显示 3更早的评论
Torsten
Torsten 2024-3-6
编辑：Torsten 2024-3-6
According to your code, your boundary conditions are
f(zeta=0) = 1
df/dzeta (zeta=0) = 0
Is this correct ?
Did you differentiate depsilon/dzeta somewhere to insert it into the differential equation ? I cannot find it in your code.
Babr
Babr 2024-3-6
Yes, the boundary conditions you said are correct.
The second term in right hand of equation has been usually overlooked in most of the studies in view of negligible impact.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-3-6
编辑：Torsten 2024-3-6
Ran in:
[t, y] = ode45(@odefun2,[0 3],[1,0]);
plot(t,y(:,1))
function dfdt = odefun2(x,y)
w = 2e16;
r0 = 30e-6;
c = 3e8;
p = 2.5;
Wp0 = p .*c ./r0;
a02 = 0.050;
phi2 = (Wp0/w)^2 * 3*a02/y(1)^4 / (4*(1+a02/(2*y(1)^2))^1.5) *...
(1 + 1/(p^2*y(1)^2) * (56+a02/y(1)^2)/(3*(1+a02/(2*y(1)^2))^0.5));
dfdt = [y(2); y(1) .^(-3) - (w .*r0 ./c) .^2 .*phi2 .*y(1)];
end
  5 个评论
显示 3更早的评论
Torsten
Torsten 2024-3-6
编辑：Torsten 2024-3-6
Ran in:
syms f(x) p c r0 a02 w
Wp0 = p*c/r0;
e = 1 - (Wp0/w)^2 / sqrt(1+a02/(2*f^2)) * (1- 1/p^2 * 3*a02/f^4 / sqrt(1+a02/(2*f^2)));
simplify(diff(e,x))
ans(x) = 
and df/dx in your code is y(2).
[t, y] = ode45(@odefun2,[0 3],[1,0]);
plot(t,y(:,1))
function dfdt = odefun2(x,y)
w = 2e16;
r0 = 30e-6;
c = 3e8;
p = 2.5;
Wp0 = p .*c ./r0;
a02 = 0.050;
phi2 = (Wp0/w)^2 * 3*a02/y(1)^4 / (4*(1+a02/(2*y(1)^2))^1.5) *...
(1 + 1/(p^2*y(1)^2) * (56+a02/y(1)^2)/(3*(1+a02/(2*y(1)^2))^0.5));
e = 1 - (Wp0/w)^2 / sqrt(1+a02/(2*y(1)^2)) * (1- 1/p^2 * 3*a02/y(1)^4 / sqrt(1+a02/(2*y(1)^2)));
sigma1 = a02/(2*y(1)^2) + 1;
dedx = 3*a02^2*c^2*y(2)/(r0^2*w^2*y(1)^7*sigma1^2) ...
-12*a02*c^2*y(2)/(r0^2*w^2*y(1)^5*sigma1)...
-a02*c^2*p^2*y(2)/(2*r0^2*w^2*y(1)^3*sigma1^1.5);
dfdt = [y(2); y(1)^(-3)-1/(2*e)*dedx*y(2)-(w*r0/c)^2*phi2*y(1)];
end

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by