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，计算积分： $\int^{1}_{-1}{e^{-x^2}dx}$

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

举例2，计算二重积分 $\int^{1}_{0}{dx}\int^{5}_{0}{sin(xy)dy}$

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

举例3，计算积分 $\int^{1}_{0}{e^{-x} + sin(x) dx}$

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

2 微分方程数值解

在这里插入图片描述

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，解方程
$\frac{2t}{y} \\y(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(0) = 1 \\ y'(0) = 0$

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

举例3，求解方程组
$f+g\\g' = -f +g\\f(0)=1\\g(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)

在这里插入图片描述
这题用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