4th order Runge Kutta Method Differential System
33 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
b = abs(27.673314419-(2*tan(1.50)))
b = 5*sin(2*pi*50*0.017)
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
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.
Antworten (1)
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
0 Kommentare
Siehe auch
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!