Solving with ode45 with a solution dependant variable in the ode

1 Ansicht (letzte 30 Tage)
Hugo
Hugo am 7 Aug. 2024
Kommentiert: Hugo am 7 Aug. 2024
Hi,
Let's say I want to solve :
with , with all following the same ode as y(z).
I can theoretically solve each ode with RK4 at step z+dz and then calculate N(z+dz). But how do I do it using the already implemented RK4 function implemented in Matlab? Would something like this work or is there a way better solution?
while(z<zLim)%%This is just an idea of the algorithm I would use
zspan=[0 dz];%Solving the ODE for a very small zspan
if(z~=0)
y0=y(length(y));
else
y0=0;
end
[z,y] = ode45(@(z,y) N*y, zspan, y0);
%Calculate c(z),g(z),....
%Calculate N at step dz
z=z+dz;
end
Thanks

Akzeptierte Antwort

Torsten
Torsten am 7 Aug. 2024
Bearbeitet: Torsten am 7 Aug. 2024
fun = @(t,y) sum(y)*y;
zspan = [0 3e-2];
y0 = [1 2 3 4 5];
[T,Y] = ode45(fun,zspan,y0);
plot(T,Y)

Weitere Antworten (1)

Aquatris
Aquatris am 7 Aug. 2024
Bearbeitet: Aquatris am 7 Aug. 2024
No, you should not do it that way. Instead you should create your function properly in a way that represents the equations. It will be easier to manipulate things later on.
[z,y] = ode45(@myfun, zspan, y0);
function dx = myfun(z,y)
% x = [y(z) N(z) g(z) c(z)]'
persistent z_prev
if isempty(z_prev)
z_prev = 0; % assuming zspan starts at 0
end
dz = z-z_prev;
dx = [x(1)*x(2); % dy = y(z)*N(z)
(x(1)+x(3)+x(4)-x(2))/dt; % dN = (y(z)+c(z)+g(z)-N)/dz so that N(z+dz) = N(z)+(y(z)+c(z)+g(z)-N(z))/dz*dz = y(z)+c(z)+g(z)
...% dc -> formulate similar to dN
...% dg -> formulate similar to dN]
end
  1 Kommentar
Hugo
Hugo am 7 Aug. 2024
I do not get how you reach this formula for dx? I tried to reach the same expression of N(z+dz) but I did not manage to do so :
Differentiating I get
Leading to :
But I can not get rid of the z+dz terms, wich are not defined at this step? The only I can see myself reach something close to your definition is by writing :
Leading to
Sorry but I am mathematically lost?
Also why not define dc and dg as dy instead of dN?

Melden Sie sich an, um zu kommentieren.

Tags

Produkte


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by