用matlab编程实现四阶龙格库塔解二元二阶微分方程组

急！！！求matlab 用四阶龙格-库塔法求解常微分方程~

建立.m文件
---------------------------------------------
function theta=danbai(t,X)
x=X(1);
dx=X(2);
ddx=-sin(x);
theta=[dx;ddx];
----------------------------------------------
命令窗口输入
>> [t,Y]=ode45(@danbai,[0 6],[pi/3 -1/2]);
>> plot(t,Y(:,1),'ro-',t,Y(:,2),'bv-');
>> legend('heta-t','dheta-t')
自编龙格库塔
------------------------------------------------
function [y,z]=Runge_kutta(a,b,y0,z0,h)
x=a:h:b;
y(1)=y0;
z(1)=z0;
n=(b-a)/h+1;
for i=2:n
K(1,1)=f1(x(i-1),y(i-1),z(i-1));
K(2,1)=f2(x(i-1),y(i-1),z(i-1));
K(1,2)=f1(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2);
K(2,2)=f2(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2);
K(1,3)=f1(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2);
K(2,3)=f2(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2);
K(1,4)=f1(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h);
K(2,4)=f2(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h);
y(i)=y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4));
z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4));
end
y(2)
plot(y,'r') %θ-t图
hold on
plot(f1(x,y,z),'g') % dθ-t图
hold on
plot(f2(x,y,z),'b-')%d2θ-t图
%f1.m
function f1=f1(x,y,z)
f1=z;
%f2.m
function f2=f2(x,y,z)
f2=-sin(y);
-----------------------------------------
Runge_kutta(0,6,pi/3,-1/2,0.02)

http://wenku.baidu.com/view/fa7ee0ebaeaad1f346933fb7.html 这里有个三阶龙格库塔的程序，可以作为参考

求解二阶微分方程，初始条件还需要给出y1'(0)和y2'(0)。这里暂时按照0处理。

 

function zd530003514

a=0.1;

b=0.1;

Y0 = [b-1; 0; b; 0];

 

% 解方程

[t,Y]= ode45(@ode,[0 10],Y0);

y1=Y(:,1);

y2=Y(:,3);

 

% 绘图 

subplot 211

plot(t,y1);

subplot 212

plot(t,y2);

 

% 微分方程定义 

function dY = ode(t, Y)

L1=5;

L2=0.01;

a0=2;

b0=2;

c0=2;

 

y1=Y(1);y2=Y(3);

dY = [

    Y(2);

    -(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1 - L1^2*y1;

    Y(4);

    -(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1;

];

 

function dy=fun(t,y)
a0=2;
b0=2;
c0=2;
lamda1=5;
lamda2=0.01;
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=-(a0*y(2)+b0*y(2)^2+c0*y(2)^3)-lamda1^2*lamda2*y(1)-lamda1^2*y(1);
dy(3)=y(4);
dy(4)=-(a0*y(2)+b0*y(2)^2+c0*y(2)^3)-lamda1^2*lamda2*y(1);
[T,Y]=ode45('fun',[0 100],[0 0 0.1 0]);

matlab解微分方程，函数调用的形式是一样的，把方程转换为代码就可以了。四阶龙格库塔应该是众多解法中的一个。详细参阅help文件。

你也弄毕设呢？

是的，所以很急
我也用MATLAB做毕设，不过你这个比我的可难多了

---

用matlab编程实现四阶龙格库塔解二元二阶微分方程组视频