how can i resolve this equation Runge kutta method

1 回表示 (過去 30 日間)
adem ski
adem ski 2020 年 1 月 22 日
編集済み: James Tursa 2020 年 1 月 23 日
y''+1=0
ddy/ddt + 1 =0
how can i resolve this equation with runge kutta method
  6 件のコメント
James Tursa
James Tursa 2020 年 1 月 22 日
@adem: You need to show us what you have done so far and then we can help you with your code.
adem ski
adem ski 2020 年 1 月 22 日
clear all
h=0.01;
t=0:0.01:10;
y(1)=1;
m=1;
f=inline('............');
for t=0:0.01:10
k1=h*feval(f,t,y(m));
k2=h*feval(f,t+(1/2)*h,y(m)+(1/2)*k1);
k3=h*feval(f,t+(1/2)*h,y(m)+(1/2)*k2);
k4=h*feval(f,t+h,y(m)+k3);
y(m+1)=y(m)+(1/6)+(k1+2*k2+2*k3+k4);
m=m+1;
end
i don't know how to define the equation?

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

回答 (2 件)

Jim Riggs
Jim Riggs 2020 年 1 月 22 日
Ths is the equation for a parabola.
y'' = -1
y' = -y
y = -1/2 y^2
No need for a numerical approximation.
  2 件のコメント
adem ski
adem ski 2020 年 1 月 22 日
i want use Runge Kutta i don't know how can i Define the equation
Jim Riggs
Jim Riggs 2020 年 1 月 22 日
ddy = @(t, y) -1;

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


James Tursa
James Tursa 2020 年 1 月 22 日
編集済み: James Tursa 2020 年 1 月 23 日
You've got a 2nd order equation, so that means you need a 2-element state vector. The two states will be y and y'. All of your code needs to be rewritten with a 2-element state vector [y;y'] instead of the single state y. I find it convenient to use a column vector for this. So everywhere in your code that you are using y, you will need to use a two element column vector instead. E.g., some changes like this:
y = zeros(2,numel(t)+1); % pre-allocate the solution, where y(:,m) is the state at time t(m)
:
y(:,1) = [1;0]; % need to initialize two elements at first time, position and velocity
:
f = @(t,y)[y(2);-1]; % [derivative of position is velocity; derivative of velocity is acceleration = -1]
:
k1 = h*feval(f, t , y(:,m) ); % changed y(m) to y(:,m)
Also, you have a bug in the line that combines the k's ... you need (1/6)* instead of (1/6)+
Make an effort at implementing these changes and then come back with any problems you continue to have.

Community Treasure Hunt

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

Start Hunting!

Translated by