Running ode45 from problem statement
3 ビュー (過去 30 日間)
古いコメントを表示
ok so my problem statement is "Use the MATLAB built-in function ode45 to perform the same calculations and generate the plots (study Example 10-11)." I have attatched the picture of the problem statement as well and I was able to run the " Sys2ODEsRK4" No problem but Im having trouble with the ode45 :
This is the Sys2ODEsRK4 code
function [t,y,ydot] = Sys2ODEsRK4(ODE1,ODE2,a,b,h,x1,y1)
t(1) = a; y(1) = .025; ydot(1) = 0;
n= (b-a)/h;
for i=1:n
t(i+1) = t(i) + h;
tm = t(i) + h/2;
Kx1 = ODE1(t(i),y(i),ydot(i));
Ky1 = ODE2(t(i),y(i),ydot(i));
Kx2 = ODE1(tm,y(i)+ Kx1*h/2,ydot(i)+Ky1*h/2);
Ky2 = ODE2(tm,y(i)+ Kx1*h/2,ydot(i)+Ky1*h/2);
Kx3 = ODE1(tm,y(i)+ Kx2*h/2,ydot(i)+Ky2*h/2);
Ky3 = ODE2(tm,y(i)+ Kx2*h/2,ydot(i)+Ky2*h/2);
Kx4 = ODE1(t(i + 1),y(i)+ Kx3*h,ydot(i)+Ky3*h);
Ky4 = ODE2(t(i + 1),y(i)+ Kx3*h,ydot(i)+Ky3*h);
y(i+1) = y(i)+(Kx1 + 2*Kx2 + 2*Kx3 + Kx4)*h/6;
ydot(i+1) = ydot(i)+(Ky1 + 2*Ky2 + 2*Ky3 + Ky4)*h/6;
end
Now for my script file part two I got the correct graphs
clc; clear all
a=0;b=10;h=0.02;
[t,y,ydot]=Sys2ODEsRK4(@dydt,@dydotdt,a,b,h,0.025,0);
figure(1)
plot(t,y,'LineWidth',2)
xlabel('Time(s)')
ylabel('Displacement (m)')
title('Displacement over time')
grid on
figure(2)
plot(t,ydot,'LineWidth',2)
xlabel('Time(s)')
ylabel('Velocity (m/s)')
title('Velocity vs Time')
grid on
function dydx=dydt(t,y,ydot)
dydx=ydot;
end
function dydx=dydotdt(t,y,ydot)
L=0.1;g=9.80;
dydx=(-2*g/L)*y;
end
Now when I do my ode45 part this is where im getting my error in the line 6 of this code:
%% Solution using matlab ode45() function
timespan = [0 10];
initCond = 0;
[t,y,ydot] = ode45(@dydt,@dydotdt,timespan,initCond);
%% plotting results
figure
plot(t,y(:,1));
xlabel('time(s)');
ylabel('Displacement (m)');
title('ode45: Displacement vs time');
figure
plot(t,y(:,2));
xlabel('time(s)');
ylabel('Velocity (m/s)');
title('ode45: Velocity vs time');
function dydx=dydt(t,y,ydot)
dydx=ydot;
end
function dydx=dydotdt(t,y,ydot)
L=0.1;g=9.80;
dydx=(-2/L*g)*y;
end
Now this is the example in the textbook for ode45 if this helps but Im lost on what im doing wrong:
function dNdt = PopRate(t,N)
bG = 1.1; bL = 0.00025; dG = 0.005;dL = 0.7;
f1 = bL*N(1)*N(2)-dL*N(1);
f2 = bG*N(2) - dG*N(2)*N(1);
dNdt = [f1;f2];
tspan = [0 25];
Nini = [ 500 3000];
[Time Pop] = ode45(@PopRate, tspan, Nini);
plot(Time,Pop(:,1),'-',Time,Pop(:,2),'--')
xlabel('Time (yr)')
ylabel('Population')
legend('Lions','Gazelles')
0 件のコメント
採用された回答
Peter O
2020 年 4 月 17 日
Hi Daniel,
For your problem, you were very close and had the solutions to each equation worked out correctly. The function signature for ode45 does not spit back the differential. This is a situation where you can take advantage of formatting the differential equations as a vector and have MATLAB solve them simultaneously. Essentially, you have a two state system, where one state is the derivative of the other, but no state needs to be integrated more than once at any timestep. That is:
dx = f(x), where x is a 2-element vector having x(1) equal to y and x(2) as y_dot. Then the solution becomes:
y = x(1);
ydot = x(2);
dx(1) = ydot;
dx(2) = -2gy/L
This is a very powerful trick that is used routinely in solving complex systems.
%% Solution using matlab ode45() function
timespan = [0 10];
initCond = [0.025, 0]; % y, ydot, grouped together
[t,y] = ode45(@dydt,timespan,initCond);
%% plotting results
figure
plot(t,y(:,1));
xlabel('time(s)');
ylabel('Displacement (m)');
title('ode45: Displacement vs time');
figure
plot(t,y(:,2));
xlabel('time(s)');
ylabel('Velocity (m/s)');
title('ode45: Velocity vs time');
function dx=dydt(t,x)
% Compute both derivatives. The states (y and ydot) come in as the x vector, of size 2x1 for time t
L=0.1;
g=9.80;
y = x(1);
ydot = x(2);
dx = zeros(2,1);
dx(1) = ydot;
dx(2) =(-2/L*g)*y;
end
4 件のコメント
Peter O
2020 年 4 月 17 日
It seems to run OK here from a single file using MATLAB's ability to nest a function within a script. Is it potentially calling an old version of dydt? The variable x is undefined in your original code. Type
which dydt
and make sure it's executing the file you expect.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Function Creation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!