【数学建模】14 微分方程模型求解方法

系统管理员 2022-12-27 14:11 308阅读 0赞

### 目录 ###

*  1 MATLAB数值微积分
 *  2 微分方程数值解
 *  3 MATLAB求解常微分方程
 *  4 课后习题

# 1 MATLAB数值微积分 #

（1）差分与微分  
• taylor 符号泰勒展开  
• polyder 多项式求导  
• diff 数值差分或符号求导

dx = diff(x) %返回向量x的差分

• gradient 数值梯度

Fx = gradient(F,x) % 返回向量F表示的一元函数沿x方向的数值梯度（F'(x)）,其中，x是与同维度的向量.x是自变量，F是因变量
    [Fx,Fy] = gradient(F,x,y)%返回矩阵F表示的二元函数的数值梯度(F'x,F'y),当F为m*n矩阵时x、y分别为n维和m维的向量。x和y是因变量

举例，求梯度，两种方式  
第一种：通过x和y的差分比.误差比较大  
第二种：用gradient函数求。两端的数据偏大，中间的数据比较准确

clear
    x = [1 1.1 1.2 1.3];
    y = x.^3;
    %第一种方式
    dy = diff(y)./diff(x) 
    %第二种方式
    dy = gradient(y,x)

（2）积分  
• trapz：梯形法积分

Z = trapz(x,y)% x是表示积分区间的离散化向量，y是与x同维数的向量，表示被积函数；z返回积分的近似值。

• dblquad：矩形区域的二重积分(二重积分)

Z = dblquad('Fun',a,b,c,d)
    %其中：Fun表示被积函数的M函数名，a、b表示x的上下限，c,d表示x的上下限

• quad：变步长数值积分

Z = quad('Fun',A,B,Tol)
    % 其中：Fun表示被积函数的M函数名，A、B上下限，Tol精度，缺省值为1e-3

• int：符号积分

int(s) %符号表达式s的不定积分
    int(s,v)% 符号表达式s关于变量v的不定积分
    int(s,a,b)%表达式s的定积分，a，b分别为下限和上限

• quad8：高精度数值积分  
举例1，计算积分： z = ∫ − 1 1 e − x 2 d x z = \\int^\{1\}\_\{-1\}\{e^\{-x^2\}dx\} z=∫−11e−x2dx

clear
    x = -1:0.1:1;
    y = exp(-x.^2);
    % 方法一
    trapz(x,y)
    % 方法二
    quad(@(x)exp(-x.^2),-1,1)

举例2，计算二重积分 z = ∫ 0 1 d x ∫ 0 5 s i n ( x y ) d y z = \\int^\{1\}\_\{0\}\{dx\}\\int^\{5\}\_\{0\}\{sin(xy)dy\} z=∫01dx∫05sin(xy)dy

Z = dblquad(@(x,y)sin(x.*y),0,1,0,5)

举例3，计算积分 z = ∫ 0 1 e − x + s i n ( x ) d x z = \\int^\{1\}\_\{0\}\{e^\{-x\} + sin(x) dx\} z=∫01e−x\+sin(x)dx

t2 = int(exp(-x) + sin(x),0,1)
    t2 = vpa(t2) %求近似解

# 2 微分方程数值解 #

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70_pic_center]

# 3 MATLAB求解常微分方程 #

• ode23 二、三阶Runge-Kutta法  
• ode15s 刚性方程组解法  
• ode45 四、五阶Runge-Kutta法  
• dsolve 符号解析解

[tout,yout] = ode45('yprime',tspan,y0)
    %这是字符串yprime是用以表示f(t,y)的M文件名
    %tspan = [t0,tf] 表示自变量初始值t0和终值tf.y0表示初值向量y0。
    %输出列向量tout表示节点（t0,t1,..,tn）
    %输出矩阵yout表示数值解，每一列对应y的一个分量。

举例1，解方程  
 y ′ = t − 2 t y y ( 0 ) = 1 , 0 < t < 1 y' = t- \\frac\{2t\}\{y\} \\\\y(0) = 1,0<t <1 y′=t−y2ty(0)=1,0<t<1

function f = fun(t,y)
    f = y-2*t./y
    f = f(:);
    end
    [t,y] = ode45('fun',[0,1],1);
    plot(t,y)

举例2，解高阶方程  
 y ′ ′ = c c o s ( 2 x ) − y y ( 0 ) = 1 y ′ ( 0 ) = 0 y'' = ccos(2x)- y \\\\ y(0) = 1 \\\\ y'(0) = 0 y′′=ccos(2x)−yy(0)=1y′(0)=0

fun = dsolve('D2y = cos(2*x)-y','y(0)=1','Dy(0)=0','x');%D2表示二阶导数，D表示一节导数
    simplify(fun)%简化的方式

举例3，求解方程组  
 f ′ = f + g g ′ = − f + g f ( 0 ) = 1 g ( 0 ) = 2 f' = f+g\\\\g' = -f +g\\\\f(0)=1\\\\g(0)=2 f′=f\+gg′=−f\+gf(0)=1g(0)=2

S = dsolve('Df=f+g','Dg =g-f','f(0)=1','g(0)=2');%求出的S是一个结构体
    S.f
    S.g

# 4 课后习题 #

计算下列微分方程组

> u’=0.5u(w-u)/v  
> v’=0.5(w-u)  
> w’=\[0.9-1000(w-y)-0.5w(w-u)\]/z  
> z’=0.5(w-u)  
> y’=-100(y-w)  
> 界条件为u(0)=v(0)=w(0)=1,z(0)=-10,w(1)=y(1)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70_pic_center 1]  
这题用dsolve去求解本人认为思路是对的，但是这样输出的结果是S =\[ empty sym \]  
有可能答案不对，如果有知道的小伙伴欢迎指出。

S = dsolve('Du=(0.5)*u*(w-u)/v','Dv =(-0.5)*(w-u)','Dw=[0.9-1000*(w-y)-(0.5)*w*(w-u)]/z',...
     'Dz=0.5*(w-u)','Dy=-100*(y-w)','u(0)=1','v(0)=1','w(0)=1','z(0)=-10','w(1)=y(1)')%求出的S是一个结构体

还有一种思路是用ode45去求精确值，但是本人不知道如何代入w(1)=y(1)这个条件。给出程序模型，表示本人的思路，如果后面明白了，会补充上。

%主程序
     x0=[1;1;1;-10;？];%初始条件
     [t,yy]=ode45(@funtest,[0,10],x0);%10是随便给出的，该括号内内容表示t的取值区间
     
     %函数
    function S=funtest(t,x)
    u=(0.5)*x(1)*(x(3)-x(1))/x(2);
    v=-0.5*(x(3)-x(1));
    w=[0.9-1000*(x(3)-x(5))-(0.5)*x(3)*(x(3)-x(1))]/x(4);
    z=0.5*(x(3)-x(1));
    y=-100*(x(5)-x(3));
    S=[u;v;w;z;y];
    end

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70_pic_center]: /images/20221120/6d2fec02311a486fade6affcf5b5e50b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkzNTY5Ng_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221120/4cd5b376790443f3a567d895083779a9.png