常微分方程欧拉法代码_欧拉法求解常微分方程

2024-08-27 11:21:03

1.如何用matlab计算微分方程的初值问题

2.关于改进欧拉法计算常微分方程，急！

3.大佬救命这matlab题怎么做？

程序已经写了，不过步长你得自己调，当步长较小时，计算时间会很长

另外，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('与理论值误差')

如何用matlab计算微分方程的初值问题

3.6.2

龙格-

库塔方法

改进的欧拉法比欧拉法精度高的原因在于，它在确定平均斜率时，多取了一个点的斜

率值。这样，如果我们在[Xi,X(i+1)]上多取几个点的斜率值，然后对它们作线性组合得到平均

斜率，则有可能构造出精度更高的计算方法。这就是龙格-库塔法的基本思想。龙格-库塔

法可看作是欧拉法思想的提高，属于精度较高的单步法。

龙格-库塔法是求解常微分方程初值问题的最重要的方法之一。MATLAB中提供了几

个用龙格-库塔法来求解常微分方程的函数，即ode23，ode45，ode113

，ode23s

，ode15s

等，其中最常用的函数是

ode23(

二三阶龙格-库塔函数)和ode45(

四五阶龙格-库塔函数)，

下面分别对它们进行介绍。

．二三阶龙格-

库塔函数(ode23)

函数

ode23

的调用格式如下：

(1)

[T,Y]=ODE23('F',TSPAN,Y0)

输入参数中的'F'

是一个字符串，表示微分方程的形

式，也可以是

)的M

文件。TSPAN=[T0

TFINAL]表示积分区间，Y0表示初始条件。

函数

ode23

表示在初始条件

Y0下从

T0到TFINAL

对微分方程

'(,)

yFty

进行积分。函数

F(T,

必须返回一列向量，两个输出参数是列向量

与矩阵

Y，其中向量

包含估计响应

的积分点，而矩阵

的行数与向量

的长度相等。向量

中的积分点不是等间距的，这是

为了保持所需的相对精度，而改变了积分算法的步长。为了获得在确定点T0,T1,

"的解，

TSPAN=[T0

TFINAL]

。需要注意的是：TSPAN中的点必须是单调递增或单调递减的。

(2)

[T,Y]=ODE23('F',TSPAN,Y0,OPTIONS)

其中，参数

options

为积分参数，它可由函

数ODESET

来设置。Options参数最常用的是相对误差‘RelTol’(

默认值是

1e-3)和绝对误差

‘AbsTol’(默认值是

1e-6)，其他参数同上。

(3)

[T,Y]=ODE23('F',TSPAN,Y0,OPTIONS,P1,P2,…)

参数P1,P2,

…可直接输入到函数

中去．如

F(T,Y,FL,P1,P2,…)。如果参数

OPTIONS为空，则输入

OPTIONS=[

]。也可

以在

ODE文件中(可参阅

ODEFILE函数)指明参数

TSPAN、Y0和OPTIONS的值。如果参

数TSPAN

或Y0

是空，则ODE23函数通过调用ODE文件[TSPAN,

Y0,

OPTIONS]

F([

],[

'init

')来获得

ODE23函数没有被提供的自变量值。如果获得的自变量表示空，则函

数ODE23会忽略，此时为

ODE23('F')。

(4)

[T,Y,TE,YE,IE]=ODE23('F',TSPAN,Y0,OPTIONS)

此时要求在参数

options

中的事

件属性设为'on'

，ODE文件必须被标记，以便

P(T,Y,'events')

能返回合适的信息，详细可参

阅函数

ODEFILE。输出参数中的

TE是一个列向量，矩阵

YE的行与列向量

TE中元素相

对应，向量

表示解的索引。

．四五阶龙格-

库塔函数(ode45)

函数

ode45

的调用格式同

ode23

相同，其差别在于内部算法不同。如果'F'

为向量函数，

则ode23

和ode45

也可用来解微分方程组。

例3.47

分别用二三阶龙格-库塔法和四五阶龙格-库塔法解常微分方程的初值问题：

解：先将微分方程写成自定义函数

exam2fun.m

function

f=exam2fun

(x,y)

f=-y-x*y.^2;

f=f(:);

然后在命令窗口输入以下语句：

[x1,y1]=ode23('exam2fun',[0:0.1:1],1)

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

0.9006

0.8046

0.7144

0.6314

0.5563

0.4892

0.4296

0.3772

0.3312

0.2910

[x2,y2]=ode45('exam2fun',[0:0.1:1],1)

0.1000

0.2000

0.3000

0.4000

0.5000

0.6000

0.7000

0.8000

0.9000

1.0000

0.9006

0.8046

0.7144

0.6315

0.5563

0.4892

0.4296

0.3772

0.3312

0.2910

关于改进欧拉法计算常微分方程，急！

我只有欧拉法可以吗。我有前欧拉法和后欧拉法。你要哪个呢？

先给你前欧拉法

function [tF,yF] = ForwardEuler(f,a,b,h,y0)

N = fix((b - a)/h);

t=a %set initial values

y=y0

tF = zeros(N+1,1); % initialize tF and yF with zero arrays

yF = zeros(N+1,length(y));

tF(1) = t; % store initial values in arrays tF, yF

yF(1,:) = y.'; % Note: y values are stored as ROW of array yF

for i=1:N

y = y + h*feval(f,t,y); % evaluate slope= feval(f,t,y);

t = t + h;

tF(i+1,:) = t; % store new values in arrays tF, yF

yF(i+1,:) = y.'; % Note: y values are stored as ROW of array yF

end

运行之后调入你的函数就ok了。a，b是范围，h是步长，y0是初值

大佬救命这matlab题怎么做？

由y'=y得y=ce^x

设y=c(x)*e^x

代入原方程

则c'(x)=(x+1)/e^x

则c(x)=-(x+1)e^(-x)-e^(-x)+c

因此，y=[-(x+1)e^(-x)-e^(-x)+c)e^x=-x-2+ce^x

把y(0)=0代入得c=2

因此，y=-x-2+2e^x

问题常微分方程数值解问题。用预估校正Euler法，求解初值问题。

求出步长h=0.1的所有点的值，并绘制图形。

求解方法用预估校正Euler法来求解，其方法是：

第一步，根据y(0)=1边界值，通过折线法计算，提供初值，即

上述式（1）也就是预报公式。

第二步，根据初值，通过梯形法计算，得到较精确的值，即

上述式（2）也就是校正公式。

这里，yn—表示y(xn)的近似值；h=x(i+1)-x(i)—表示步长

第三步，按上述循环计算，计算当x分别等于0，0.1，0.2，0.3，0.4，0.5，0.6，0.7，0.8，0.9，1的y(x)值。

第四步，根据x和y值，进行绘制该微分方程的数值解曲线。

求解过程

通过计算，得到微分方程的数值解

x =[0 0.1 0.20.30.4 0.50.6 0.7 0.8 0.91.0]

y =[1 1.0959 1.1841 1.26621.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379]

根据这些点，就可以绘制其图形。

matlab求解求解步骤：

第一步，根据预估校正Euler法的迭代式，编写其函数，如Euler_Cauchy（func,x0,y0,xf,h）

这里，func——表示微分方程，x0,y0——表示微分方程初始值，xf——表示x的终值，hf——表示步长

第二步，编写微分方程函数，即

function f=func(x,y);?

f(1)=y-2*x/y;

f=f(:);

end

第三步，编写主程序，即

y0=1; ?%初值

x0=0;xf=1; %x的范围

n=10; %等份

h=(xf-x0)/n %步长

[t,x]=Euler_Cauchy(@func,x0,y0,xf,h); %预估校正Euler法计算

plot(x,y,'ks-') %绘制函数图

tu on %绘制坐标区网格线

xlabel('x'),ylabel('y(x)') %标注坐标轴名称

第四步，运行后，即可得到如下图形。

数值解图形分析从上述图形中，可以看到得到的数值解与精确解比较，精度较高的顺序依次是：

欧拉法→改进的欧拉法（预估校正Euler法）→龙格－库塔法→解析解

常微分方程数值解方法

1、欧拉法。欧拉方法（也叫折线法）是最早的一种数值方法。欧拉方法是一种数值解微分方程的方法，它是由瑞士数学家欧拉发明的。欧拉方法的基本思想是将微分方程转化为差分方程然后通过迭代求解差分方程来逼近微分方程的解。是一种一阶数值方法，用以对给定初值的常微分方程（即初值问题）求解。它是一种解决数值常微分方程的最基本的一类显型方法。

欧拉方法的具体步骤如下：

首先将微分方程转化为差分方程，即将微分方程史的导数用差分代替，然后将差分方程史的未知函数值用前一时刻的函数值代替，得到一介递推公式。接着.从初始时刻开始，按照递推公式依次计算出每个时刻的函数值，直到达到所需的时刻为止。最后，将计算出的函数值作为微分方程的近似解。

欧拉法迭代式

2、预估校正欧拉法。预估校正欧拉法是对欧拉算法的改进方法。微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。

这个方法中，(1)式用折线法提供初值，称为预报公式。(2)式用梯形法给出较精确的值，称为校正公式。合称预报校正公式。