发布网友
共2个回答
热心网友
程序已经写了,不过步长你得自己调,当步长较小时,计算时间会很长
另外,tend是时间的终值,你可以设小一些。因为解析解为10*cos(x),我设成pi,就是计算半个周期。
x''(t)=-x(t)
引入y1=x,y2=x',则
y1'=y2
y2'=-x=-y1
初始条件为:
y1(0)=10;
y2(0)=0;
将下面两行百分号之间的内容,保存成DiffEulerRk4.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MaxDiffX=DiffEulerRk4(dt,PlotFlag)
%dt是步长
%PlotFlag是否作图
if nargin<1
dt=0.01;
end
if nargin<2
PlotFlag=0;
end
f=inline('[y(2);-y(1)]','t','y'); %微分方程的右边项
t0=0; %初始时刻
tend=pi; %计算的点数
tt=t0:dt:tend; %一系列离散的点
N=length(tt); %点数
y0=[10;0];
%%(1)欧拉法
EulerY=y0;
for i=2:N
EulerY(:,i)=EulerY(:,i-1)+dt*f(tt(i-1),EulerY(:,i-1));
end
EulerX=EulerY(1,:); %取x
%%(2)龙格库塔法
RkY=y0;
for i=2:N
k1=f(tt(i-1), RkY(:,i-1));
k2=f(tt(i-1)+dt/2, RkY(:,i-1)+k1*dt/2);
k3=f(tt(i-1)+dt/2, RkY(:,i-1)+k2*dt/2);
k4=f(tt(i-1)+dt, RkY(:,i-1)+k3*dt);
RkY(:,i)=RkY(:,i-1)+(k1+2*k2+2*k3+k4)*dt/6;
end
RkX=RkY(1,:); %取x
%精确解
syms t
analytic=dsolve('D2x=-x','x(0)=10','Dx(0)=0','t');
rightdata=subs(analytic,t,tt);
if PlotFlag
plot(tt,EulerX,'b-',tt,RkX,'r--',tt,rightdata,'g-.')
legend('Euler','Runge-Kutta','analytic')
end
MaxDiffX=[max(abs(RkX-rightdata)),max(abs(EulerX-rightdata))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
所有题,都得你自己调步长。
输入:
DiffEulerRk4(0.01,1) %步长取0.01的计算结果,参数为1代表作图,自己得修改步长
%%下面是变化
Error=[];
Dt=[5e-4,1e-3,2e-3,5e-3,0.01,0.05,0.1];
for dt=Dt %几种步长,自行修改
dt %查看dt,步长小,计算量大
Error_1=DiffEulerRk4(dt); %不作图
Error=[Error;Error_1]; %保存欧拉法误差
end
semilogx(Dt,Error)
legend('Euler','RK4')
xlabel('步长')
ylabel('误差')
title('与理论值误差')
请参考
热心网友
要建立项目的,2010你直接新建一个文件是不可以的,这个可能也是很多人用VC++ 6.0的原因,直接新建文件就可以编译了。VC2010先要新建项目,再在项目里添加文件