Question

大神你好 我又来问你matlab的问题啦 请赐教

上次你给我解答的问题在这里：
http://zhidao.baidu.com/question/523737127?quesup2&oldq=1
这个问题里面现在等于已经把x1，x2，u，m，lambda1，lambda2求出来了。
但是我还想实现一个功能，及时求积分，积分的形式在图中。（我自己试过什么int函数quad函数都不行，请大神帮忙啊）

Answer 1

1、这个积分不能用int或quad之类的函数来求。使用这些函数的前提是，你要有被积函数关于积分变量t的解析表达式，而事实上，被积函数只有通过数值方法计算出来离散点，没有表达式，所以不能用int或quad函数求。

其实这个积分可以直接根据概念来计算：现在已经有了函数 f 在区间上一些离散点 tk 处的函数值 fk（k=1..n），那么，函数 f 对变量 t 在该区间的积分大致可以表示为（矩形法近似）
Σfk*Δtk
具体到函数代码，就是在上次问题的追问部分计算出u和m之后，加上下面几句：
x1 = sol.y(1,:);x2 = sol.y(2,:);
f=50*(x1-24).^2+60*(x2-0.75).^2+60*(u-40).^2+30*(m-50).^2;
J=1/2*sum(f(1:end-1).*diff(t))
这样计算得到的结果是 4.6119e+004。

2、上面只是根据概念可以想出的最简单做法，。为了提高精度，可以考虑采用以下措施：
（1）使用梯形法求解。
计算数值积分，最容易想到的就是使用梯形法代替矩形法，也就是把上面三行代码的最后一行改成
J=1/2*sum((f(1:end-1)+f(2:end))/2.*diff(t))
这样计算出来的结果是 4.5712e+004。

（2）使用插值。
很容易想到，数值积分的精度受Δt影响很大，越小就越接近 dt 的概念，也就越精确（可以用 min(diff(t)) 和 max(diff(t)) 查看BVP求解结果的Δt范围）。这样，可以考虑使用更小的Δt，例如
T=0:0.001:4;
F=interp1(t,f,T);
J=1/2*sum((F(1:end-1)+F(2:end))/2.*diff(T))
矩形法和梯形法得到的结果比较接近，分别为 4.5713e+004 和 4.5712e+004。注意，这里的插值默认使用分段线性插值，对提高精度帮助不大，如果使用样条或三次插值，会得到更好的结果（均为4.5733e+004）。

（3）BVP解的初始估计把点取密集一些。
插值属于求解BVP之后的后处理方法，另一种思路是，求解BVP的时候就把点取密集一些。例如把原代码在最前面（解的初始估计）改为
solinit = bvpinit(linspace(0,4,1000),[5 0.5 1 2]);
使用梯形法求解得到结果是 4.5733e+004。

3、我不知道你所说的“及时求积分”是什么意思？因为积分上限是4，而且就问题本身而言属于边值问题，所以必须在完整求解微分方程之后才可能计算该积分，而不可能在解微分方程的过程中（比如，在t=1的时刻）计算。

4、请再仔细核对一下表达式是否正确。我回答的重点在方法，具体表达式不保证录入会不会有错。

花了一个多小时的时间，希望对你有帮助。