Question

微分方程数值解法，知道导数的边界值，而不是函数的边界值，怎么做？

我知道微分方程边值问题，把微分变成差分，然后解方程组可解。但是这个问题边界值不是函数的边界值，而是函数导数的边界值，请问怎么处理？看了2本书，也没答案

比如y''+y'+y=f(x) a1<=x<=a2
y'(a1)=b1,y'(a2)=b2

Answer 1

这个方程最简单的方法是用shooting method.可wiki 之。
主要就假设一个y(a1), 然后看，y'(a2)是否等于b2。
可以用matlab写两个小函数，一个用在ode45（fun1），然后y(a1)的值可以用fzero(fun2).
fun1就是把二阶微分化成一阶方程组，fun2就是调用ode45并返回y'(a2)-b2. 这样会自动求出初值的。

别外这个方程可能有解析解。可以利用格林函数来解，得到一个积分形式的解。

给你看一个我以前写的例子吧。适当改下就行了。
function out=fdif(eta,f)
out=zeros(size(f));
out(1)=f(2);
out(2)=f(3);
out(3)=f(2)^2-1-f(1)*f(3);

function result=soldif(x)
[etaout,fout]=ode45('fdif',[0,10],[0,0,x]);
result=fout(end,2)-1;

追问

您好，最近看了一下shooting method，发现我这个问题似乎不能应用此办法，我的方程齐次项是y''+9/4y=0，shooting method是把原方程分解为y=u+c*v，c是常数，v就是齐次方程。x的区间是[-pi,pi]，这就导致了一个问题，因为shooting method要求v‘（-pi）=0，所以v’（pi）=0，c无法求解。我看的shooting method例子边界条件是y=。。。而不是y‘=。。。请问我该如何处理，两边求导？

追答

你对shooting method的理解好像不太对(或者说和我说的不一样)。对于齐次方程，比如你说的y''+9/4y=f(x)。应该像我之前说的先变成一次微分方程组。y1'=y2; (y2)'+9/4y1=f(x)所以下面函数改成
function out=fdif(eta,f)
out=zeros(size(f));
out(1)=f(2);
out(2)=f(x)-9/4*f(1); //注意这个f(x)改成相应的，如f(x)=x^2; then 改成eta^2
这样就可以调用ode45.
把下面代码改成
function result=soldif(x)
[etaout,fout]=ode45('fdif',[-pi,pi],[x,0]); // [x,0]分别为假设初如值及其一阶导数
result=fout(end,2)-0; // 这里的0是最后一点导数。
最后利用fzero就可以求解，试试吧。Good luck

追问

您好，解二阶初值问题微分方程的方法我知道，也就是ode45之前我都知道。但是我看书上和wiki百科，关于线性问题，y=y1(t)+(y1-y1(t))/2*y2(t)，（篇幅太长，您可以自己上wiki上看看）。但是我这个问题y2(t)=0，所以我就不知道如何处理了。wiki百科上第一部分说的和您一样，是要解方程的根，但是我不知道是解什么方程的跟，他上边写的是:F(a)=y(t1;a)-y1,其中y(t0)=y0,y(t1)=y1,假设初值是y(t0)=y0,y'(t0)=a。

追答

wiki那个问题是已经知道两个端点的值。所以先假设一个端点的值y0,一阶导数值为a,这样求解程后，会得到另一个端点的值，和已知的端点值y1比较，直到差为0.我上面替你改的程序可以直接用的。只要改下f(x)。为什么这样我大体说了。fzero就是求那个F(a)的。

追问

您还是没明白我的意思，我取了个特殊情况，能求出解析解，发现解里确实有个不确定常数，所以shooting method里的(y1-y1(t))/2*y2(t)确实求不出来，换成常数就能和解析解对应上。THX，anyway

Answer 2

这个是线性非其次方程，解出来的通解一定是 齐次式通解加特解形式。
特解，可以用你的数值解法来找一个。
这样就能得到原方程的通解。通解中有两个待定常数。
把y'(a1)=b1,y'(a2)=b2代入进去，得到两个方程。可以把两个常数解出来

追问

大哥，您搞笑吗，那个要是有解析解，谁求数值解。我的问题是怎么求数值解，谢谢

追答

y''+y'+y=0这主用一般方法可以解的吧
y''+y'+y=f(x),用数值解法可以找特解的吧。(比如幂级数解法)

凑在一起，解可以找出来数值解的。我晕。

还有一个办法 两边求导
得到y'''+y''+y' =f'(x)

自己没明白别人的意思，还说明人高笑，不想细讲了。

追问

恩。两边求导可以，但是太麻烦了，不过是机器算，也无所谓了

追答

第一个办法是本质。好好领会，不要学了半天方程，就背了几个破公式。记得采纳。

追问

请问幂级数法不是求解其次非线性方程的吗，怎么求解特解，有没有其他方法求解特解？推荐本书看看，谢谢

追答

谁告诉你幂级数集法只能求齐次方程的。
只要方程满足足够的解析条件，都能解的。
你在任何一本微分方程的教材，都有相关的介绍。

Answer 3

使用matlab用龙格-库塔解法求数值解，可以调用函数，也可以利用simulink做。
简单的用matlab的符号运算可以得到解析解。

有限元中一般使用伽辽金法（余量加权法的一种），简单的例子可参考http://19se.com.cn/a/gongju/2011/0511/515.html
解微分方程组对工程问题来说十分重要，正不断学习中。

追问

您好，龙格库塔是求解初值问题的，而不是边值问题，ode45也是，谢谢