求matlab 用四阶龙格-库塔法求解常微分方程
来源:学生作业帮 编辑:搜狗做题网作业帮 分类:综合作业 时间:2024/07/22 18:25:54
求matlab 用四阶龙格-库塔法求解常微分方程
用matlab~
2θ/dt2 +sin θ=0 0
用matlab~
2θ/dt2 +sin θ=0 0
![求matlab 用四阶龙格-库塔法求解常微分方程](/uploads/image/z/17580025-1-5.jpg?t=%E6%B1%82matlab+%E7%94%A8%E5%9B%9B%E9%98%B6%E9%BE%99%E6%A0%BC-%E5%BA%93%E5%A1%94%E6%B3%95%E6%B1%82%E8%A7%A3%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B)
建立.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('\theta-t','d\theta-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://img.wesiedu.com/upload/6/fd/6fd8ce42625aa693f7cf97b8ab9bb75e.jpg)
---------------------------------------------
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('\theta-t','d\theta-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://img.wesiedu.com/upload/6/fd/6fd8ce42625aa693f7cf97b8ab9bb75e.jpg)