二阶偏微分方程组 龙格库塔法_常微分法数值计算方法

fb1bd624b1b60d5600168539eee8f56d.png

假设一阶常微分方程组有下式给出

其中矢量

为状态变量
构成的向量,即
,常称为系统的状态向量,n称为系统的阶次,而
为任意函数数,t为时间变量,这样就可以采用数值方法求解常微分方程组了。另外任意高阶微分方程都可以通过变量替换变成一阶微分方程组,这里不再赘述.

求解常微分方程组的数值方法是多种多样的,如常用的Euler法、Runge-Kutta方法、Adams线性多步法、Gear法等。

  1. 欧拉法
    时刻系统状态向量的值为
    ,若选择计算步长h,则可以写出在
    时刻系统状态向量的值为

    这样,用迭代的方法可以由给定的初值问题逐步求出在所选择的时间段t∈[0,T]内各个时刻
    处的原问题数值解。提高数值解精度的一种显然的方法是减小步长h的值。然而,并不能无限制地减小h的值,这主要有两条原因:

    (1)减慢计算速度.因为对选定的求解时间而言,减小步长就意味着增加在这个时间段内的计算点数目,故计算速度减慢;
    (2)增加累积误差因为不论选择多小的步长,所得出的数值解都将有个舍入误差,减小计算步长则将增加计算的次数,从而使得整个计算过程的舍入误差的叠加和传递次数增多,产生较大的累积误差。
  2. 四阶龙格库塔法(Runge-Kutta)

3. Adams算法

matlab常微分方程求解器

ode求解器​ww2.mathworks.cn

求解非刚性微分方程 ,ode45,ode23,ode15s,ode223等等,用啊基本类似,以ode45为例说明.

51cd12a897fe00619f3ec28a041da436.png

语法

[t,y] = ode45(odefun,tspan,y0)
[t,y] = ode45(odefun,tspan,y0,options)
[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
sol = ode45(___)

说明示例
[t,y] = ode45(odefun,tspan,y0)(其中 tspan = [t0 tf])求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分,初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

所有 MATLAB® ODE 求解器都可以解算 y′=f(t,y) 形式的方程组,或涉及质量矩阵 M(t,y)y′=f(t,y) 的问题。求解器都使用类似的语法。ode23s 求解器只能解算质量矩阵为常量的问题。ode15s 和 ode23t 可以解算具有奇异质量矩阵的问题,称为微分代数方程 (DAE)。使用 odeset 的 Mass 选项指定质量矩阵。

ode45 是一个通用型 ODE 求解器,是您解算大多数问题时的首选。但是,对于刚性问题或需要较高准确性的问题,其他 ODE 求解器可能更适合。有关详细信息,请参阅选择 ODE 求解器。

[t,y] = ode45(odefun,tspan,y0,options) 还使用由 options(使用 odeset 函数创建的参数)定义的积分设置。例如,使用 AbsTol 和 RelTol 选项指定绝对误差容限和相对误差容限,或者使用 Mass 选项提供质量矩阵。

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) 还求 (t,y) 的函数(称为事件函数)在何处为零。在输出中,te 是事件的时间,ye 是事件发生时的解,ie 是触发的事件的索引。

对于每个事件函数,应指定积分是否在零点处终止以及过零方向是否重要。为此,请将 'Events' 属性设置为函数(例如 myEventFcn 或 @myEventFcn),并创建一个对应的函数:[value,isterminal,direction] = myEventFcn(t,y)。有关详细信息,请参阅 ODE 事件位置。

sol = ode45(___) 返回一个结构体,您可以将该结构体与 deval 结合使用来计算区间 [t0 tf] 中任意点位置的解。您可以使用上述语法中的任何输入参数组合。

举例

[t,y] = ode45(@vdp1,[0 20],[2; 0]);
% Plot the solutions for $y_1$ and $y_2$ against |t|.
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')
function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];


c3beac4ed933df6c82b356082a54eb9d.png
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 1024 设计师:白松林 返回首页