Solving with ode45 with a solution dependant variable in the ode

3 vues (au cours des 30 derniers jours)
Hugo
Hugo le 7 Août 2024
Commenté : Hugo le 7 Août 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

Réponse acceptée

Torsten
Torsten le 7 Août 2024
Modifié(e) : Torsten le 7 Août 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)

Plus de réponses (1)

Aquatris
Aquatris le 7 Août 2024
Modifié(e) : Aquatris le 7 Août 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 commentaire
Hugo
Hugo le 7 Août 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?

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Tags

Produits


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by