差商近似1阶导数matlab,常微分方程的解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似...

常微分方程的解法求解系列博文:

常微分方程的解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似

常微分方程的解法 (二): 欧拉(Euler)方法

常微分方程的解法 (三): 龙格—库塔(Runge—Kutta)方法 、线性多步法

常微分方程的解法 (四): Matlab 解法

目录

1 常微分方程的离散化

数值解法

(i)用差商近似导数------差分方程初值问题

(ii)用数值积分方法

(iii)Taylor 多项式近似

建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并 加以检验。如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线 性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得到这样的解, 而绝大多数变系数方程、非线性方程都是所谓“解不出来”的,即使看起来非常简单的 方程如 9d0bf6df61f57e96b1415f5d78b5e6a1.png ,于是对于用微分方程解决实际问题来说,数值解法就是一个十 分重要的手段.

1 常微分方程的离散化

下面主要讨论一阶常微分方程的初值问题,其一般形式是

a5b4705ac5ef667182438cdec70ed228.png

在下面的讨论中,我们总假定函数 f (x, y) 连续,且关于 y 满足李普希兹(Lipschitz)条 件,即存在常数 L ,使得

352f4683030c4a07a25d3f9589218e14.png

这样,由常微分方程理论知,初值问题(1)的解必定存在唯一。

数值解法

所谓数值解法,就是求问题(1)的解 y(x) 在若干点

a15a11333d9886e0f5c5376d0d263ece.png

建立数值解法,首先要将微分方程离散化,一般采用以下几种方法:

(i)用差商近似导数------差分方程初值问题

691ac3f83f78aeaafc1a7dcfa92ac817.png

需要说明的是,用不同的差商近似导数,将得到不同的计算公式。

(ii)用数值积分方法

将问题(1)的解表成积分形式,用数值积分方法离散化。例如,对微分方程两端 积分,得

18ee50bac893014f11a5a9608994199e.png

右边的积分用矩形公式或梯形公式计算。

(iii)Taylor 多项式近似

e1b091554192ef40d3e3ecac8a3c75c2.png

以上三种方法都是将微分方程离散化的常用方法,每一类方法又可导出不同形式的 计算公式。其中的 Taylor 展开法,不仅可以得到求数值解的公式,而且容易估计截断 误差。

常微分方程的解法求解系列博文:

常微分方程的解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似

常微分方程的解法 (二): 欧拉(Euler)方法

常微分方程的解法 (三): 龙格—库塔(Runge—Kutta)方法 、线性多步法

常微分方程的解法 (四): Matlab 解法

  • 1
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
matlab程序 matlab s函数 . function [sys,x0,str,ts] = sfuntmpl(t,x,u,flag) % SFUNTMPL 是M-文件 S函数模板 % 通过剪裁,用户可以生成自己的S函数,不过一定要重新命名 % 利用S函数可以生成连续、离散混合系统等,实现任何模块的功能 % % M-文件S函数的语法为: % [SYS,X0,STR,TS] = SFUNC(T,X,U,FLAG,P1,...,Pn) % % 参数含义: % t是当前时间 % x是S函数相应的状态向量 % u是模块的输入 % flag是所要执行的任务 % % FLAG 结果 功能 % ----- ------ -------------------------------------------- % 0 [SIZES,X0,STR,TS] 模块初始 % 1 DX 计算模块导数 % 2 DS 更新模块离散状态 % 3 Y 计算模块输出 % 4 TNEXT 计算下一个采样时间点 % 9 [] 结束仿真 % % % 用户切勿改动输出参数的顺序、名称和数目 % 输入参数的数目不能小于1,这四个参数的名称和排列顺序不能改动 % 用户可以根据自己的要求添加额外的参数,位置依次为第5,6,7,8,9等。 % S函数的flag参数是一个标记变量,具有6个不同值,分别为0,1,2,3,4,9 % flag的6个值分别指向6个不同的子函数 % flag所指向的子函数也成为回调方法(Callback Methods) switch flag, %初始,调用“模块初始”子程序% case 0, [sys,x0,str,ts]=mdlInitializeSizes; %连续状态变量计算,调用“计算模块导数”子函数% case 1, sys=mdlDerivatives(t,x,u); %更新,调用“更新模块离散状态”子函数% case 2, sys=mdlUpdate(t,x,u); %输出,调用“计算模块输出”子函数% case 3, sys=mdlOutputs(t,x,u); %计算下一时刻采样点,调用“计算下一个采样时刻点”子函数% case 4, sys=mdlGetTimeOfNextVarHit(t,x,u); %结束,调用“结束仿真”子函数% case 9, sys=mdlTerminate(t,x,u); %其他的flag% otherwise DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag)); end % end sfuntmpl %============================================================================= % “模块初始”子函数 % 返回大小、初始条件和样本 function [sys,x0,str,ts]=mdlInitializeSizes % 调用simsizes函数,返回规范格式的sizes构架 % 这条指令不要修改 sizes = simsizes; % 模块的连续状态个数,0是默认值 % 用户可以根据自己的要求进行修改 sizes.NumContStates = 0; % 模块的离散状态个数,0是默认值 % 用户可以根据自己的要求进行修改 sizes.NumDiscStates = 0; % 模块的输出个数,0是默认值 % 用户可以根据自己的要求进行修改 sizes.NumOutputs = 1; % 模块的输入个数,0是默认值 % 用户可以根据自己的要求进行修改 sizes.NumInputs = 2; % 模块中包含的直通前向馈路个数,1是默认值 % 用户可以根据自己的要求进行修改 sizes.DirFeedthrough = 1; % 模块中采样时间的个数,1是默认值,至少需要一个样本时间 % 用户可以根据自己的要求进行修改 sizes.NumSampleTimes = 1; % 初始后的构架sizes经simsizes函数处理后向sys赋值 % 这条指令不要修改 sys = simsizes(sizes); % 给模块初始值变量赋值,[]是默认值 % 用户可以根据自己的要求进行修改 x0 = []; % 系统保留变量 % 切勿改动,保持为空 str = []; % “二元对”描述采样时间及偏移量。[0 0]是默认值 % [0 0]适用于连续系统 % [-1 0]则表示该模块采样时间继承其前的模块采样时间设置 ts = [0 0]; % end mdlInitializeSizes %============================================================================= % 计算导数向量 function sys=mdlDerivatives(t,x,u) % 此处填写计算导数向量的指令 % []是模块的默认值 % 用户必须把算得的离散状态的导数向量赋给sys sys = []; % end mdlDerivatives %============================================================================= % 计算离散状态向量 function sys=mdlUpdate(t,x,u) % 此处填写计算离散状态向量的指令 % []是模块的默认值 % 用户必须把算得的离散状态向量赋给sys sys = []; % end mdlUpdate %============================================================================= % 计算模块输出向量 function sys=mdlOutputs(t,x,u) % 此处填写计算模块输出向量的指令 % []是模块的默认值 % 用户必须把算得的模块输出向量赋给sys global m l J g=9.8; M0=1/3*m*l*l; G0=1/2*m*l*g*cos(u(1)); % 期望轨迹 dd_qr=-sin(t);d_qr=cos(t);qr=sin(t); d_q=u(2);q=u(1); % u1=(M0+J)*(dd_qr+2*2*(d_qr-d_q)+4*(qr-q))+G0; sys = u1; % end mdlOutputs %============================================================================= % 计算下一采样时刻 function sys=mdlGetTimeOfNextVarHit(t,x,u) % 该子函数仅在“采样时间”情况下使用 % sampleTime = 1 是模块默认设置,表示在当前时间1秒后再调用本模块 % 用户可以根据自己的要求修改 sampleTime = 1; % 将计算得到的下一采样时刻赋给sys % 切勿改动 sys = t + sampleTime; % end mdlGetTimeOfNextVarHit %============================================================================= function sys=mdlTerminate(t,x,u) % 模板默认设置,一般情况不要改动 sys = []; % end mdlTerminate

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值