-2

eaxmple first-order derivative using matlab.if we did for second-order derivative,How does this code change?

 % y'=(x+2)/y    y(0)=1
 % y=sqrt(x^2+4*x+1)
 clear all
 x(1)=0;

 ry(1)=1;
 dx=0.1;
  for i=1:100
k1=(x(i)+2)/ry(i);
k2=(x(i)+dx/2+2)/(ry(i)+dx/2*k1);
k3=(x(i)+dx/2+2)/(ry(i)+dx/2*k2);
k4=(x(i)+dx+2)/(ry(i)+dx*k3);
ry(i+1)=ry(i)+dx*(k1+2*k2+2*k3+k4)/6;

 x(i+1)=x(i)+dx;
 end
 y=sqrt(x.^2+4*x+1);
 plot(x,y,'r-',x,ey,'b-',x,hy,'y-',x,ry,'g-');
user3713186
  • 1
  • 1
  • Try to put your differential equation inside a function, the code may become much clearer. Make it generic enough that it can work with vectors in place of the x(i). – Lutz Lehmann Jun 06 '14 at 12:25

1 Answers1

0

maybe this code help you if this is your topic and problem, it is runge kutta code for solving second order derivative problems,this belong to 2 years ago and i wrote it by some help.this problem is defined in the book i said in the last answer. 

syms x y1 y2       
f=input('input function like:D2y=(-0.1*y2)-x, y1 & y2 subtitute as y & Dy respectively      and do not type D2y in front of f input');    

f=inline(f);    
g=inline(y2);           
a1=0;   
b1=1;
x1=0;  
h=0.5;  
n=4;  
A=zeros(n,3);  
A(1,1)=0;  
A(1,2)=a1;  
A(1,3)=b1;   


for i=1:n      
%%%%%%%%%%%%%%k1    
F1 = h*(subs(g,[x,y1,y2],[x1,a1,b1]));   
F2 = h*(subs(f,[x,y1,y2],[x1,a1,b1]));   
W1=F1;   
W2=F2;    
%%%%%%%%%%%%%%%%k2   
a2=a1 + (F1/2);    
b2=b1 + (F2/2);    
x2= x1 + h/2;     
F1 = h*(subs(g,[x,y1,y2],[x2,a2,b2]));     
F2 = h*(subs(f,[x,y1,y2],[x2,a2,b2]));      
W1= W1+ (2*F1);     
W2= W2+ (2*F2);    
%%%%%%%%%%%%%%%k3    
a3=a1 + (F1/2);     
b3=b1 + (F2/2);    
x3=x1 + h/2;     
F1 = h*(subs(g,[x,y1,y2],[x3,a3,b3]));     
F2 = h*(subs(f,[x,y1,y2],[x3,a3,b3]));      
W1= W1+ (2*F1);     
W2= W2+ (2*F2);     
%%%%%%%%%%%%$$k4     
a4=a1 + (F1);     
b4=b1 + (F2);     
x4=x1 + h;    
F1 =h*(subs(g,[x,y1,y2],[x4,a4,b4]));     
F2 = h*(subs(f,[x,y1,y2],[x4,a4,b4]));     
W1= W1+ (F1);    
W2= W2+ (F2);     
%%%%%%%%%%%%%%%%    
W1 = (W1/6);    
W2 = (W2/6);     
a1=a1+ W1;     
b1=b1+ W2;    
x1=x1 + h;    
A(i+1,1)=x1;    
A(i+1,2)=a1;    
A(i+1,3)=b1;    
end       
fprintf('y(2) of runge kutta: %15.14f\n',A(n+1,2));    
fprintf('Dy(2)of runge kutta: %15.14f\n',A(n+1,3));
m.j
  • 1