matlab实现欧拉法和RK-4方法的数值计算

Question

已知系统模型和初始条件如下: x''(t)=-x(t),x(0)=10,x'(0)=0

（1）采用C语言等高级语言编制程序，步长取0.01，积分方法采用欧拉法和四阶龙格-库塔。改变仿真步长，观察仿真结果，找出系统由稳定到不稳定的临界步长。
（2）你认为该系统步长取多大时精度最高。
（3） 利用MATLAB仿真软件，求解上述线性定常系统的数值解，要求输出y随步长变化的值（最好能够以图形显示）。

要有全部程序,悬赏直接200,有的发邮件到[email protected],编译通过立即给分

Answer 1

程序已经写了，不过步长你得自己调，当步长较小时，计算时间会很长
另外，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('与理论值误差')

Answer 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=[];
title('与理论值误差'

Answer 3

　　一、计算方法：
　　unction [x,y]=DEEuler(f,a,b,y0,n);
　　%f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end
　　h=(b-a)/n;
　　x(1)=a;y(1)=y0; for i=1:n
　　x(i+1)=x(i)+h;
　　y(i+1)=y(i)+h*feval(f,x(i),y(i)); end
　　format short

　　二、MATLAB的简单介绍：
　　MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
　　MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

Answer 4

这题好难啊！