如何在Matlab中用四阶龙格库塔法解微分方程
3个回答
展开全部
建立.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)
---------------------------------------------
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)
展开全部
function
dy=test(t,y)
dy=[-12*cos(y(2))-120*cos(208*2*pi/360-3*y(2));(12*sin(y(2))+120*sin(208*2*pi/360-3*y(2)))/y(1);];
[t,y]=ode45('test',[0.01,1],[1,1])
plot(t,y(:,1),t,y(:,2));
x=3000*sin(70*2*pi/360)-y(1).*sin(y(2));
z=3000*cos(70*2*pi/360)-12*t-y(1).*cos(y(2));
plot(t,x,t,z)
最后的这图像不好。你可再看看。方法是没问题的。
你可到matlab中文论坛去请教关于matlab的问题,那儿高手多的很。
dy=test(t,y)
dy=[-12*cos(y(2))-120*cos(208*2*pi/360-3*y(2));(12*sin(y(2))+120*sin(208*2*pi/360-3*y(2)))/y(1);];
[t,y]=ode45('test',[0.01,1],[1,1])
plot(t,y(:,1),t,y(:,2));
x=3000*sin(70*2*pi/360)-y(1).*sin(y(2));
z=3000*cos(70*2*pi/360)-12*t-y(1).*cos(y(2));
plot(t,x,t,z)
最后的这图像不好。你可再看看。方法是没问题的。
你可到matlab中文论坛去请教关于matlab的问题,那儿高手多的很。
已赞过
已踩过<
评论
收起
你对这个回答的评价是?
展开全部
dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
y=(3*sin(5*x))/exp(2*x)
一楼有错误
y=(3*sin(5*x))/exp(2*x)
一楼有错误
已赞过
已踩过<
评论
收起
你对这个回答的评价是?
推荐律师服务:
若未解决您的问题,请您详细描述您的问题,通过百度律临进行免费专业咨询