Solving with ode45 with a solution dependant variable in the ode

2 ビュー (過去 30 日間)
Hugo
Hugo 2024 年 8 月 7 日
コメント済み: Hugo 2024 年 8 月 7 日
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

採用された回答

Torsten
Torsten 2024 年 8 月 7 日
編集済み: Torsten 2024 年 8 月 7 日
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)
  1 件のコメント
Hugo
Hugo 2024 年 8 月 7 日
Thanks that is way smarter than my strategy

サインインしてコメントする。

その他の回答 (1 件)

Aquatris
Aquatris 2024 年 8 月 7 日
編集済み: Aquatris 2024 年 8 月 7 日
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 件のコメント
Hugo
Hugo 2024 年 8 月 7 日
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?

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

製品


リリース

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by