Filter löschen
Filter löschen

4th order Runge Kutta Method Differential System

33 Ansichten (letzte 30 Tage)
Prince Nino
Prince Nino am 17 Sep. 2024 um 19:51
Bearbeitet: Torsten am 17 Sep. 2024 um 23:38
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 Kommentare
Torsten
Torsten am 17 Sep. 2024 um 19:56
Bearbeitet: Torsten am 17 Sep. 2024 um 19:59
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 am 17 Sep. 2024 um 19:57
dont mind
"a = 0.025/2500
b = abs(27.673314419-(2*tan(1.50)))
b = 5*sin(2*pi*50*0.017)"

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Torsten
Torsten am 17 Sep. 2024 um 23:27
Bearbeitet: Torsten am 17 Sep. 2024 um 23:38
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

Kategorien

Mehr zu Manage Products finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by