常偏微分方程解高手请参与讨论!

bzzz

matlab的ode45里关于变步长部分写的有点复杂，它考虑了很多情况
下面这个思想差不多，但是思路清晰很多

Max=10;
Min=1e9;
% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/Max;
hmin = (tfinal - t)/Min;
h = (tfinal - t)/100; % initial guess at a step size
......%some statements are skipped
% estimate the local truncation error
gamma1 = x5 - x4;
% Estimate the error and the acceptable error
delta = norm(gamma1,'inf');       % actual error
tau = tol*max(norm(x,'inf'),1.0); % allowable error

% Update the solution only if the error is acceptable
if (delta<=tau),
         t = t + h;
         x = x5;   
% <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
end % if (delta<=tau)

      % Update the step size
if (delta==0.0),
       delta = 1e-16;
end % if (delta==0.0)
h = min(hmax, 0.8*h*(tau/delta)^pow);
%one method for variable step calculation

可以看到
变步长未必会需要更长的时间
因为它是一种自适应的方式
如果某阶段情况理想的话
就会用最大步长来做
那就会很快
按初始的step guess要100步
如果变步长情况良好，最快只要10步^^

BTW，上面只是一种变步长的方式，可以用很多不同的，关键是整个流程和思路
简言之，就是
用误差范数和atol来看精度，进而来确定下一个步长

hehe
如有疏漏，请指正^_^

damingsun

lyrock  你好！

看到你的解答，很受启发，在此先谢了。
我现在学matlab主要用来求解声学的一些问题，与你的advanced vibration很像。以后还要多向你请教呢。呵呵

bzzz

今天在研学上看到的一个问题

[parse]heigher wrote:[/parse]
模型：图1
单质点弹簧系统，质量M，弹簧刚度K，施加的力为F，初始位移、速度均为0。求解系统在F作用下动力响应，主要是得到质量块位移随时间变化曲线。
在MATLAB中编写如下程序求解上述微分方程
function odesolve(n)
clear;
clc;
M=10;
K=500;
F=8;
tspan=[0:1:700]; ％求解区间t
[t,y]=ode23(@sol,tspan,[0;0],[ ],M,K,F);  ％采用ode23求解
plot(t,y(:,1),'-')
%----------------------------------------------
function dydt = sol(t,y,M,K,F)         ％微分方程
dydt=[y(2);-K/M*y(1)+F/M];
得到的结果：图2

问题：因为不考虑任何摩擦和阻尼，求解该问题的结果应当是一正弦曲线的样子，不会随时间衰减的，而现在衰减了，请教问题出在什么地方？另外图片怎样加上去？我怎么加不上去

首先说明一个问题
这样写tspan=[0:1:700];
只会让ode45/23作定步长求解，在你所指定的700个点求解
造成波形的误差
然后对于后面的波形问题
如果用tspan=[0,200],ode45会用变步长
但是做出的图形就会像你那样，出现包络线是一条斜线
这是由于数值解的误差的传播
解决方法是减小误差容限
options = odeset(RelTol',10e-7)
[t,y1]=ode45(@sol,[0,200],[0;0],options,M,K,F);
可以看到变化
不过tspan取得这么长没什么意义，这么简单的方程不如用符号函数
下图蓝色的是减小误差容限，红色的是默认容限

朱顶红

请教高手们：我在用ODE113解一非线性微分方程组时，得到的数值解是复数。用PLOT绘图时只用实部。出现这样的警告：Warning: Imaginary parts of complex X and/or Y arguments ignored.我深思不得其解。请指点！谢谢！

hncscp306

热和流体领域中需要处理的基本上都是偏微分方程,真正能用解析法求解的少的可怜,只能是利用数值方法求解.最主要的数值方法有两种,有限差分法和有限单元法.

damingsun

runge-kutta法是针对初值问题的，对于边值问题，大家一般处理？是不是象我在上面的帖子里说的，要定义一个方程：Dy（2）＝1。
突然明白了一下，现在的理解是：如果是边值问题，直接把x当成t处理好了，这主要是对我们这些新手来说的，对于老手可能根本就不存在二者的区别。
请高手评论！

bainhome

FEMLAB版有个置顶的贴子,有关PDE的,转一下,大家可以看看:
http://www.simwe.com/forum/viewthread.php?tid=240745
另外关于楼上所问关于边值问题,本来早想和大家探讨,无奈这两天太忙,我先给一个老薛的例子:
求解线性方程边值问题的打靶法的例子
  function [t,y]=shooting(f1,f2,tspan,x0f,varargin)
  t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2);
  [t,y1]=ode45(f1,tspan,[1;0],varargin);
  [t,y2]=ode45(f1,tspan,[0;1],varargin);
  [t,yp]=ode45(f2,tspan,[0;0],varargin);
  m=(gb-ga*y1(end,1)-yp(end,1))/y2(end,1);
  [t,y]=ode45(f2,tspan,[ga;m],varargin);
示例:
求解线性微分方程D2y-3*Dy+2*y=x,y(0)=1,y(1)=2在[0,1]段方程的数值解.
解:两个函数
      function xdot=c7fun1(t,x)
      xdot=[x(2); -2*x(1)+3*x(2)];
     function xdot=c7fun2(t,x)
     xdot=[x(2); t-2*x(1)+3*x(2)];
求解:
[t,y]=shooting('c7fun1','c7fun2',[0,1],[1;2]);
绘制解图:
plot(t,y)

bainhome

昨天在mathworks上看到的,界面做得十分标准,大家可以看看.这是截图:

bainhome

主程序:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6703&objectType=file
另外一个求解类似"f(t,y,y')=0 where y'=dy/dt"的程序,针对R10,估计有点老:
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=28&objectType=file

yycgy

提个小建议，希望对大家有帮助。

算法的东西，在Fortran和C中可以找到很多的源程序，如果大家有兴趣，不妨可以把这些程序用matlab改写过来或者通过接口直接调用。

个人观点，仅供参考。

bainhome

yycgy的看法从工科实用角度上是成立的,而现实中我在做数值计算程序时也奉行这种拿来式的思想,一句话:省时省力.
但这个主题有所不同:我的想法从一开始就不是为了省力气而做的,只是想把MATLAB作为"数学"软件,成为一个载体来学习一些微分方程的基本原理,无论是别人做过的,或者自己原创的,汇集到这个贴子里,使得大家集思广益,当朋友们遇到类似的问题,或者正在研究生阶段学习数值分析,都可以想到在simweMATLAB版有这样一个贴子能有万一的借鉴之处,而混编根据我的看法,对于初学者来说,其难度未必或者远远大于自己理解原理编写一个定(变)步长的RK法,所以,熟悉原理是这个贴子我的初衷,很高兴能有lyrock和bzzz两位高手的参与,这个贴子还远没有完,我想随着大家的学问日深,还会有更好更多的内容充实到这个贴子里,而yycgy提到的看法也很好,也欢迎你或者混编的高手能够做这个工作并把你们的心得写到这个主题下.

yycgy

提供一个抛物线PDE组的小程序。

bainhome

以下这段代码是xue_dy在研学上发的,从代码风格上看他很有可能是一位我非常敬仰的高手,他的书我阅读数遍始终进益良多,转一下这段代码,其中eval命令的用法我很喜欢^_^
=========================================================================

这是我写的一段代码
function [tout,yout] = rk_4(odefile,tspan,y0)
t0=tspan(1); th=tspan(2);
if length(tspan)<=3, h=tspan(3);
else, h=tspan(2)-tspan(1); th=tspan(2); end
tout=[t0:h:th]'; yout=[];
for t=tout'
k1=h*eval([odefile '(t,y0)']);
k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']);
k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']);
k4=h*eval([odefile '(t+h,y0+k3)']);
y0=y0+(k1+2*k2+2*k3+k4)/6;
yout=[yout; y0'];
end
不过还是建议别用定步长的方法

kingslee

能用解析法的方程很少,有几种方法
1.把偏微分方程做离散化处理,最常用的是Galerkin离散法,取有效的截断精度,这样偏微分方程就离散为一个常微分方程组,再用常微分理论解,其中涉及矩阵相关理论
2.摄动法的理论得到了广泛应用,所有就可以在Galerkin离散后对常微分方程组摄动,解出精度很高的近似解
3.近年来,对偏微分方程直接摄动的研究更是热点,直接摄动要比离散后摄动的解更令人满意.所以可以对偏微分方程直接摄动,本人在做这方面的可以

junziyang

我觉得Matlab的几个ODE solver功能就够强大的了。ODE45是4阶5次Runge-Kutta，ODE113是多步法....这些solver自动选择步长，非常稳定，速度也非常快。

yycgy

kingslee wrote:
能用解析法的方程很少,有几种方法
1.把偏微分方程做离散化处理,最常用的是Galerkin离散法,取有效的截断精度,这样偏微分方程就离散为一个常微分方程组,再用常微分理论解,其中涉及矩阵相关理论
2.摄动法的理论得到了广泛应用,所有就可以在Galerkin离散后对常微分方程组摄动,解出精度很高的近似解
3.近年来,对偏微分方程直接摄动的研究更是热点,直接摄动要比离散后摄动的解更令人满意.所以可以对偏微分方程直接摄动,本人在做这方面的可以

能不能提供一些摄动的matlab程序，谢谢。

ccaa

怎样在ode45中把其他中间计算变量给调出来！！

解算中方程组右端
function DYDT=dydt(t,y)
f=f(y1,y2,yn);%---------f是y1,y2,yn的函数
DYDT=[...;
       ...];

y1,y2...yn可以通过ode45求出，但变量f(y1,y2,...,yn)怎么求呢？
我也用过for 循环解决这类问题，在每步循环中可以调出变量f的数值。但解算量太大，太慢！！现在我想要在保证速度的同时又调出这些变量，怎么办呢？请高人指教！！！

smtmobly

支持lyrock ：
工具都是有的！大家都是应用嘛！比如你们说的什么算法，我是学计算数学的！
我从没有去做过这些！也没有一个老师在上课的时候提的那么清楚过！
因为计算算法关注的是计算的简单方便！既然大家都是应用的
，关心的是算法之间的区别！比如你为什么用欧拉迭代，又 为什么用梯形公式！
区别在哪！语言上描述就行了！
就是精度上的差别啊
你可以用高斯公式那一定用高斯了！为什么有的能用有的不能用，
很简单啊高斯公式对光滑性要求太高啊！
所以现在有限元都用Nystorm算法啊！
当然本人一直都没有用matlab 进行计算，都是用c语言！今天来看看
高手们写的东西有点摸不清所以提这个建议！
谢谢大家

ilxy

我自己的课题之中需要解偏微分方程组，我觉得应该自己编程序来解决，自己的程序你可以了解每一个计算的细节。也可以采用不同的计算方法，比如有限差分法或者有限单元法加上迦廖金法。编程的过程是辛苦的，但是经过这样艰苦的认真的思考以后，收获也是切实的。况且随着研究的深入，会出现新的偏微分方程形式，可能会超出现成的软件的使用范围，这个时候更需要自己的主观能动性，自己编程序去解决，自己编程可以很清楚问题的算法，而我认为算法是一个问题的灵魂。对于一个偏微分方程组，确定合理的差分格式和求解次序，这个思考的过程是最重要的。那种把微分方程直接输入到某个软件中，然后就希望得到计算结果的人，我觉得实在亵渎数学的灵魂。可能最后我也会使用某个软件，但是前提是我清楚每一步计算的细节。我不会相信没有仔细思考的结果。

Iveco

原帖由 jojobi 于 2005-7-26 11:37 发表
我个人觉得惯于'偏微分方程的数值解法'回帖少的原因主要matlab的PDE工具太简单了，对于计算流体力学、计算传热学，都要求解偏微分方程组，PDE工具就无法使用，需要自己编程了。再有PDE工具只是支持三角形单元，过 ...

       我也有和这位兄弟同样的感受，在求解偏微分方程的时候，很多时候都得自己编程实现，就像楼上的仁兄所说的，自己编程可以锻炼自己的思维过程，达到真正意义上的会解方程，其实偏微分方程的数值解法主要有两种，一、差分方法；二、有限元法。
最近学了差分法，觉得先将差分格式化为矩阵形式，然后求解矩阵方程，写了个显式差分格式计算双曲线型方程的程序。特做抛砖引玉，还请大家发表看法。

function [u,u3,value,U]=xianshi(h,t,a）
% 显式差分格式计算双曲线型方程
format short g;    % 设置输出数据的精度
s=a*t/h;         % 步长
for i=1:1/h+1      % 0<=x<=1时的边界条件的数值解(共1/h+1个点)
     x1(i)=0+(i-1)*h;
     u(1,i)=0;
end
for i = 2:1/h
    u(2,i) = t*x1(i)+t^2/2*F(x1(i),0);   % 第二行的数值
end
for j = 2:1/t+1      % 0<t<=1时的边界条件的数值解(共1/t个点)
    x2(j-1) = 0+(j-1)*t;
    u(j,1) = 0;
    u(j,1/h+1) = sin(x2(j-1));
end
e = ones(1/h-1,1);
A = spdiags([e*s^2 e*2*(1-s^2) e*s^2], -1:1, 1/h-1, 1/h-1);   %三对角系数矩阵 A 的生成  
for j = 1:1/t-1
    f(1,1) = t^2*F(x1(2),x2(j+1))+s^2*u(j+1,1);
    for k = 2:1/h-2
        f(k,1)=t^2*F(x1(k+1),x2(j+1));
    end
    f(1/h-1,1)=t^2*F(x1(1/h),x2(j+1))+s^2*u(j+1,1/h+1);
    u0=u(j+1,2:1/h);
    u1=u(j,2:1/h);
    u2=A*u0'-u1'+f;
    for k = 2:1/h
        u(j+2,k) = u2(k-1,1);
    end
end
k=0.1*1/t+1;   
t1=0.1:0.1:1;
value=G(0.5,t1)';         %取出特定位置的精确解            
w=k:k-1:1/t+1;
u3=u(w,1/(2*h)+1);        % 取出与精确解对应的数值解      
U=abs(u3-value);          %两者的绝对误差
end

function fun=F(x,t)   %函数项
fun=(t^2-x^2)*sin(x*t);
end

function  val=G(x,t)   %精确解
val=sin(x*t);
end