# how to generate a code for solving 2st order differential equations using improved runge kutta 4th order method

조회 수: 2(최근 30일)
rayhan hindi 2021년 6월 18일
편집: James Tursa 2021년 6월 18일
I know who to solve 1st order differential equations using this code:
function [tp,yp] = rk4(dydt,tspan,y0,h)
if nargin < 4, error('at least 4 input arguments required'), end
if any(diff(tspan)<=0), error('tspan not ascending order'), end
n = diff(tspan)/h;
ti = tspan(1); tf = tspan(2);
tp = ti:h:tf;
yp = zeros(ti,length(tspan));
for n = 1:length(tp)
k1 = dydt(tp(n),y0);
k2 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k1*h);
k3 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k2*h);
k4 = dydt(tp(n) + h, y0 + k3*h);
yp(n,:) = y0 + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h)
ti = ti + h;
y0 = yp(n);
end
plot(tp,yp)
end
but i dont how to solve it if the differential equation was from a secend order like : y''+2y=8x(9-x)
what changes should I do to the code so it can solve this equation.
(I didnt write the code above I got help from th website)

댓글을 달려면 로그인하십시오.

### 답변(1개)

James Tursa 2021년 6월 18일
편집: James Tursa 2021년 6월 18일
The posted code actually has a bug. Also it is almost set up properly for vector solutions, but not quite. Make these changes:
yp = zeros(numel(tspan),numel(y0)); % changed (original line using ti was a bug)
yp(1,:) = y0; % added
for n = 1:length(tp)-1 % changed
k1 = dydt(tp(n) , y0 );
k2 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k1*h);
k3 = dydt(tp(n) + (1/2)*h, y0 + (1/2)*k2*h);
k4 = dydt(tp(n) + h, y0 + k3*h);
yp(n+1,:) = y0 + ((1/6)*(k1 + 2*k2 + 2*k3 + k4)*h); % changed
% ti = ti + h; % deleted ... a useless line
y0 = yp(n+1,:); % changed
end
Then it is just a matter of defining your derivative function. For your case it is 2nd order so you will need a function that returns a two element vector. Using this equation:
y''+2y=8x(9-x)
Define a two element state vector y such that:
y(1) = the original y
y(2) = y'
Then the derivatives are simply:
d y(1) /dx = dy/dx = y' = y(2)
d y(2) /dx = dy'/dx = y'' = 8x(9-x)-2y = 8*x*(9-x)-2*y(1)
Then you simply write this function coding these derivatives:
function dy = myderivative(x,y)
dy = [y(2),8*x*(9-x)-2*y(1)];
end
Where the 2nd element is obtained by solving the y''+2y=8x(9-x) equation for y'' and using y(1) for the original y.
One caveat to this is that your posted code is set up to store the state vectors as rows of the final solution, so that is how i wrote the derivative function ... i.e., it returns a row vector. However, if you use ode45( ) your derivative function would have to return a column vector instead. Something to be aware of.

댓글을 달려면 로그인하십시오.

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by