4th order Runge Kutta Method Differential System

1 次查看（过去 30 天）
Prince Nino
Prince Nino 2024-9-17
编辑： Torsten 2024-9-17
Ran in:
4th order Runge Kutta Method Differential System. I don't know what's wrong with my code and I have been trying to figure out what went wrong. I'm trying to solve a system with 3 First-Order Differential Equations.
clear all
close all
clc
COne = 1000;
CTwo = 1000;
R = 50;
L = 0.1;
Vt = 0.026;
Is = 1*10^(-8);
a = 0.025/2500
a = 1.0000e-05
b = abs(27.673314419-(2*tan(1.50)))
b = 0.5295
b = 5*sin(2*pi*50*0.017)
b = -4.0451
xprime_func = @(t,x,y,z) ((1/COne)*y-(1/(COne*R))*x);
yprime_func = @(t,x,y,z) ((1/L)*z-(1/L)*x);
zprime_func = @(t,x,y,z) ((1/CTwo)*Is*(exp((abs(10*sin(2*pi*50*t))-z)/(2*Vt))-1)-(1/CTwo)*y);
% Define time interval and step size
tmax=0.025; steps=2500; h=tmax/steps;
% Initial conditions:
x(1)=0; y(1)=0; z(1)=0; t(1)=0;
% Estimate of derivatives and marching in time.
for i=1:steps
t(i+1)=i*h;
K(1)=h*xprime_func(t(i),x(i),y(i),z(i));
L(1)=h*yprime_func(t(i),x(i),y(i),z(i));
M(1)=h*zprime_func(t(i),x(i),y(i),z(i));
K(2)=h*xprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));
L(2)=h*yprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));
M(2)=h*zprime_func(t(i)+h/2,x(i)+1/2*K(1),y(i)+1/2*L(1),z(i)+1/2*M(1));
K(3)=h*xprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));
L(3)=h*yprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));
M(3)=h*zprime_func(t(i)+h/2,x(i)+1/2*K(2),y(i)+1/2*L(2),z(i)+1/2*M(2));
K(4)=h*xprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));
L(4)=h*yprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));
M(4)=h*zprime_func(t(i)+h,x(i)+1*K(3),y(i)+1*L(3),z(i)+1*M(3));
x(i+1)=x(i)+1/6*(K(1)+2*K(2)+2*K(3)+K(4));
y(i+1)=y(i)+1/6*(L(1)+2*L(2)+2*L(3)+L(4));
z(i+1)=z(i)+1/6*(M(1)+2*M(2)+2*M(3)+M(4));
end
plot(t,x,t,y,t,z);
  2 个评论
Torsten
Torsten 2024-9-17
编辑：Torsten 2024-9-17
The Table 7.7 and figure 7.9 is what should my values for my x,y,z and it is insanely way far off.
Please include a mathematical description of your problem with equations and constants used, Table 7.7 and figure 7.9.
I would advice to solve the problem first with a sophisticated MATLAB integrator like ode45 to get a reference solution.
Prince Nino
Prince Nino 2024-9-17
dont mind
"a = 0.025/2500
b = abs(27.673314419-(2*tan(1.50)))
b = 5*sin(2*pi*50*0.017)"

请先登录，再进行评论。

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

回答（1 个）

Torsten
Torsten 2024-9-17
编辑：Torsten 2024-9-17
Ran in:
Note that COne and CTwo are given in muF - thus they should be prescribed as 1000e-6 in your code, I guess.
COne = 1000e-6;
CTwo = 1000e-6;
R = 50;
L = 0.1;
Vt = 0.026;
Is = 1e-8;
xprime_func = @(t,x,y,z) ((1/COne)*y-(1/(COne*R))*x);
yprime_func = @(t,x,y,z) ((1/L)*z-(1/L)*x);
zprime_func = @(t,x,y,z) ((1/CTwo)*Is*(exp((abs(10*sin(2*pi*50*t))-z)/(2*Vt))-1)-(1/CTwo)*y);
f = @(t,x,y,z)[xprime_func(t,x,y,z);yprime_func(t,x,y,z);zprime_func(t,x,y,z)];
F = @(t,u)f(t,u(1),u(2),u(3));
% Define time interval and step size
tmax=0.025; steps=2500; tspan=linspace(0,tmax,steps);
% Initial conditions:
u0 = [0;0;0];
[T,U] = ode15s(F,tspan,u0);
figure(1)
plot(T,U(:,2))
grid on
figure(2)
plot(T,[U(:,1),U(:,3)])
grid on

类别

在 Help Center 和 File Exchange 中查找有关 Manage Products 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by