常微分,偏微分方程的数值解法讨论 zz

redstar · 发表于 2005-8-21 14:48

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

转载自simwe 
 
最近我在整理simweMATLAB板块的旧贴时发现在用"方程"这个关键词搜索出的贴子质量低得惊人!单贴率(这里指的是问题的未解决率)几乎可以达到
70%以上,你们在这个部分看到的贴子已经是经过我过滤的了,可是还是很没有质量,没有太多,或者是基本没有让人看完感觉有所收获的内容,更别提什么相见
恨晚,醍醐灌顶的感觉了,大部分是回答者请发问者自己去看simulink的帮助,要不就help ode45,help
pdetool...感觉技术含量不是太高.因此有一个想法:如果有在偏微分和常微分的数值解法方面有一定研究的朋友们请来这里随便谈谈您在这个方面的一
些看法和见解,未必囿于软件,一个经典常用算法的理解,一个专业方向中某问题中用到的微分方程,一个收敛精度的探讨...都可以让它丰富成一个精练的贴
子,成为后来者可资借鉴的内容.我知道,借用alexqxp版主的话说:MATLAB版神佛满天,太多高手在此默默潜水汲取灵感给养并辅助完成自身的修炼
与提高,其中不乏在数学,物理方面的真正高手,请各位百忙中抽些时间,帮助后进,bainhome在这里替许多辗转苦学不得其法的朋友们说声谢谢,当然我
也知道,有关这部分的内容在牵扯算法基本原理的地方已经比较深了,数理方程,数值逼近,差分,有限元,泛函变分...没有能一个贴子说清楚的,甚至每一个
都是引无数英雄折腰或者曾经折腰,正在折腰的领域,我自己就有一些偏微分数值解法的书,但没有一本能让我轻松读上一页而没有困惑的!而单就偏微分方程来说
已经很有意思了:有些是方程的数值算法促使了问题的解决,有些是实际问题催化了新算法的产生,剪不断,理还乱,个人看法是:如果先给几个简单或者有迹可循
的例子为出发点,大家有时间的话慢慢扩展,也许能有些新的触动,以上只是自己一厢情愿的瞎想,高手勿笑. 我自己先在这里提个比较基础的问题给个思考的起始点: 变步长RUNGE-KUTTA法的MATLAB实现(算例说明) 本
贴还是那个老规矩,凡是诸如"顶","支持","好贴","对新手太有帮助啦"这样的垃圾灌水贴一律删除.对提出并解决问题的朋友给予1~3分的积分奖
励,还是老话:分值小意思,表达一个尊重原创的态度.我也很清楚,这样的方式对真正的高手不会有太大吸引力,甚至有反效果,但是请将负面看法只责怪我本人
的世俗唐突,而并不因为反感我的做事方式而影响我所谈到的问题的解决,cwit一篇文章里有句话我比较赞同:"没有师傅，就寻求一些高手的帮助。找不到高
手，就寻找朋友。没有朋友感兴趣，就培养这个群体。土壤厚重了，我们这些生物才能生长得更茂盛一些！" 就用这个作为结尾吧^_^

redstar · 发表于 2005-8-21 14:48

正餐第一道:一种特殊的二阶RUNGE-KUTTA法,但思想是来自于改进的euler法,即:提供一个y(xn+p)的预测值再计算该点斜率值K=f(x_n_plus_b,y_n_plus_b),程序如下: function fun_ode_runge_kutta %建立二阶runge-kutta公式,即:变形的euler公式: xn=0:.1:1; y0=1; y_n=[]; for i=1:length(xn)-1 K1=y0-2*xn(i)/y0; K2=y0+.5*diff(xn([1:2]))*K1-2*(xn(i)+.5*diff(xn([1:2])))/(y0+.5*diff(xn([1:2]))*K1); y_n=[y_n;y0+diff(xn([1:2]))*K2]; y0=y_n(end); end %精确的解析解 dy=dsolve('Dy=y-2*x/y','y(0)=1','x'); y_exact=subs(dy,'x',{.1:.1:1}); %精度比较 y_compare=[y_n';y_exact]; y_compare_minor=[y_compare(2,:)-y_compare(1,:)]; var_y=var(y_compare_minor'); 和上面的程序实际是一个思路,最后做一个方差比较,结果如下: y_compare= 1.0955 1.1833 1.2651 1.3419 1.4145 1.4836 1.5497 1.6131 1.6741 1.733 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 var_y=9.7261e-008 显
然可以看出,虽然只是二阶RUNGE-KUTTA法的程序,具有二阶代数精度,在本例中的精度已经相当不错了,比上贴中的改进euler法的精度1e-6
的量级已经有了改善,二阶runge-kutta公式如下,只是在这里我们选择了特殊的系数使之成为变形的euler公式而已: 
<a href="http://www.simwe.com/forum/upload/2005/07/28/58872829.jpg" target="_blank" ><img src="http://www.simwe.com/forum/upload/2005/07/28/58872829.jpg" title="Click to view full alexqxp.jpg (347 X 210)" align="middle" border="0" height="210" width="347"></a>

redstar · 发表于 2005-8-21 14:48

正餐第二道:三阶RUNGE-KUTTA法的程序 function fun_ode_runge_kutta_3 %建立三阶runge-kutta公式,使用其中的一个特例 xn=0:.1:1; y0=1; y_n=[]; for i=1:length(xn)-1 K1=y0-2*xn(i)/y0; K2=y0+.5*diff(xn([1:2]))*K1-2*(xn(i)+.5*diff(xn([1:2])))/(y0+.5*diff(xn([1:2]))*K1); K3=y0-diff(xn([1:2]))*K1+2*diff(xn([1:2]))*K2-2*(diff(xn([1:2]))+... xn(i))/(y0-diff(xn([1:2]))*K1+2*diff(xn([1:2]))*K2); y_n=[y_n;y0+1/6*diff(xn([1:2]))*(K1+4*K2+K3)]; y0=y_n(end); end %精确的解析解 dy=dsolve('Dy=y-2*x/y','y(0)=1','x'); y_exact=subs(dy,'x',{.1:.1:1}); %精度比较 y_compare=[y_n';y_exact]; y_compare_minor=[y_compare(2,:)-y_compare(1,:)]; var_y=var(y_compare_minor'); =================================================================================== y_compare= 1.0954 1.1832 1.2649 1.3416 1.4142 1.4833 1.5492 1.6125 1.6734 1.7321 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 var_y=2.1742e-010 精度再上一个层次(三阶代数精度)
 
<a href="http://www.simwe.com/forum/upload/2005/07/28/18768317.jpg" target="_blank" ><img src="http://www.simwe.com/forum/upload/2005/07/28/18768317.jpg" title="Click to view full alexqxp.jpg (582 X 309)" align="middle" border="0" height="309" width="582"></a>

redstar · 发表于 2005-8-21 14:49

I would like to make some summaries in the view of an engineer for this topic. 1. My Background First,
let me introduce my background. I am Ph.D. candidate in mechanical
engineering, and my research focuses on structural dynamics and signal
collection & processing. So, my review can not be as theoretical as
a person in mathematic dose. I will emphasize on the practical
engineering application (not on mathematics). Also, this review is
based on my understanding of this topic, so if there is something wrong
or uncompleted, hope you guys can understand and make it completed. 2. Main contents in the summary There techniques (methods) about solving differential equations in engineering will be talked about. There are: 1) Finite Difference Numerical Computation 2) Runge-Kutta Method 3) ODE function in MATLAB In my point of view, I think they are similiar. ============================================== for
the reseaon of short-of-time, I hope i can finish the summary in 3
times. Presently, only the first part is finished. Hope it could be
useful to you guys. 
 <img src="http://www.simwe.com/forum/images_zh_CN/attachment.gif" align="middle" border="0"><a href="http://www.simwe.com/forum/user/download/544622/Numerical%20Solution%20for%20Differential%20Equation.pdf" target="_blank" >Numerical Solution for Differential Equation.pdf</a> (55.31k)

redstar · 发表于 2005-8-21 14:50

A corresponding GUI program for the "Finite Difference Numerical Computation" (FDNC) is attached.
 
 <img src="http://www.simwe.com/forum/images_zh_CN/attachment.gif" align="middle" border="0"><a href="http://www.simwe.com/forum/user/download/544626/FDNC.rar" target="_blank" >FDNC.rar</a> (5.33k) 
 

I simplized the procedure for fear someone get confused and cannot keep going. <img src="http://www.simwe.com/forum/images/smiles/smile_big.gif" alt="Big Smile" width="15">
 
<a href="http://www.simwe.com/forum/upload/2005/07/29/10297567.jpg" target="_blank" ><img src="http://www.simwe.com/forum/upload/2005/07/29/10297567.jpg" title="Click to view full Taylor.JPG (658 X 285)" align="middle" border="0" height="277" width="640"></a>

redstar · 发表于 2005-8-21 14:50

the second part of my summary
is about Runge-Kutta Method (RKM). The PDF file is updated. And hope
you guys can give some feedback. thanks. A GUI program of RKM will be presented later.
 
 <img src="http://www.simwe.com/forum/images_zh_CN/attachment.gif" align="middle" border="0"><a href="http://www.simwe.com/forum/user/download/545694/merical%20Solution%20for%20Differential%20Equation%28II%29.pdf" target="_blank" >merical Solution for Differential Equation(II).pdf</a> (94.27k)

 
 

the attached file is the Runge-Kutta Method GUI program. feedbacks are welcome, thanks
 
 <img src="http://www.simwe.com/forum/images_zh_CN/attachment.gif" align="middle" border="0"><a href="http://www.simwe.com/forum/user/download/545695/Runge_Kutta%20GUI.rar" target="_blank" >Runge_Kutta GUI.rar</a> (6.2k)

redstar · 发表于 2005-8-21 14:51

正餐第三道:四阶runge-kutta法程序 function fun_ode_runge_kutta_4 %建立三阶runge-kutta公式,使用其中的一个特例 xn=0:.1:1; y0=1; y_n=[]; for i=1:length(xn)-1 K1=y0-2*xn(i)/y0; K2=y0+.5*diff(xn([1:2]))*K1-2*(xn(i)+.5*diff(xn([1:2])))/(y0+.5*diff(xn([1:2]))*K1); K3=y0+.5*diff(xn([1:2]))*K2-2*(xn(i)+.5*diff(xn([1:2])))/(y0+.5*diff(xn([1:2]))*K2); K4=y0+diff(xn([1:2]))*K3-(2*xn(i)+2*diff(xn([1:2])))/(y0+diff(xn([1:2]))*K3); y_n=[y_n;y0+1/6*diff(xn([1:2]))*(K1+2*K2+2*K3+K4)]; y0=y_n(end); end %精确的解析解 dy=dsolve('Dy=y-2*x/y','y(0)=1','x'); y_exact=subs(dy,'x',{.1:.1:1}); % 精度比较 y_compare=[y_n';y_exact]; y_compare_minor=[y_compare(2,:)-y_compare(1,:)]; var_y=var(y_compare_minor'); ================================================================================= var_y=2.8864e-012 精
度再往上走两个量级,runge-kutta法是最常用的单步高精度微分方程的解法,ode45的基本思想即来自于此,由于lyrock对这个方法的基本
问题已经总结的比较全面,因此,我在这里只是简单介绍一下我自己学习的感受,同样,未必全面,未必精确,甚至未必正确,只是作为前面我所参与内容的一点心
得,一家之言而已: runge-kutta法的基本内涵是来自于函数的taylor展开的,我们知道taylor展开是函数逼近的有效方法.但由
于展开过程中求导的机时代价过高,因此采用了在一个步长范围内再做内插节点序列,求数值差分得到内插节点序列的差分值,然后将这些斜率平均之,由此来形成
runge-kutta法的迭代格式,显然从基本思路中可以看出,runge-kutta法这样一种高精度的单步方法仅仅可用于微分方程的初值问题,不能
解决边值问题,而所谓的二阶,三阶,四阶的runge-kutta法只是在等步长中内插入不同数量的序列节点而已,我们知道,内插点越多,平均值就越接近
真实步长经过点的导数值,只是我们在实际操作过程中多了一个机时和精确度两者之间的权衡而已.而正是因为这个权衡引出了另外一个问题:如何能保证在一定机
时的情况下精确度的最理想化?即转化为处理两个问题: 1.怎样衡量和检验计算结果的精度? 2.如何依据所获得的精度处理步长? 这样引出了变步长的runge-kutta法,我并不清楚ode45是怎样的一种变步长的runge-kutta法,但数值分析书中对这个问题的描述是: 先从节点xn出发,以h为步长求得近似值,记做y_n+1_(h),由于公式局部阶段误差为o(h^5)则有: y(x_n+1)-y_n+1_(h)≈ch^5 再将步长折半,跨两步到x_n+1求得近似值y_n+1_(h/2),每跨一步截断误差为o((h/2)^5),所以有 y(x_n+1)-y_n+1_(h)≈2c(h/2)^5 这样步长折半误差可以减少到1/16,对于最后我们所设定的精度ε,如果折半两次计算结果的偏差Δ>ε,则再折半,直到Δ<ε为止 反之,如果如果折半两次计算结果的偏差Δ<ε,则加倍步长,直到Δ>ε为止,这时再将步长折半一次,得到所要的结果. 这
样做的目的无非是得到合理的步长值使得计算速度在给定精度下最为合理,因此我们在ode45中事前并不能知道(除非你很有经验,那是另外一个故事了
^_^)步长确切等于多少(我个人的理解是这样,不知道正确与否,请大家给予指正)由于变步长的原因,我上面的工作的目的也就显现出来了,即定步长只能够
使精度达到1e-12的量级,而同样的问题,用ode45变步长方法则能够得到1e-18的精度量级,所以总体上讲虽然计算量增加,但是是合算的.其他的
问题我也等待lyrock给出更加全面的解答,他的工作应该还有最后一个部分,很期待...^_^我最近在做论文的前期编程工作,很忙,只能这样凭借记
忆,写些很粗糙的东西,上面不对的地方还请lyrock和bzzz见谅...

redstar · 发表于 2005-8-21 14:51

Damingsun　你好，大家好： 首先，MATLAB的功能实在太强大，我仍然是个学习者。学习经验谈不上，就谈点对学习MATLAB的看法吧。 对
这一观点“matlab水平不高的人搞计算，用的是工具箱。高手搞计算，自己写程序”我持否定观点。因为MATLAB本身就是来简化计算的，有方便的不
用，为什么要用复杂的？前日看到一贴子，认为“高手写GUI程序用代码直接写，入门的人写GUI程序用GUI工具箱”，对这一点我也不同意，GUI工具箱
能使你更快地实现你的目的，为什么就入门的人用？我不相信自己写代码比用GUI工具箱快。 但大家可能会问：为什么MATLAB有ODE函
数不用，还写什么RUNGE－KUTTA程序？首先，我想解释一下，这是偶上一门“Advanced
Vibration”的课时写的程序。其次，自己写程序有个好处在于：你能真正理解“what’s going on in
it”；如果在你的研究中，需要你用C语言或其它语言写这种算法，那你就会轻车熟路。有一点应该相信：学会的东西在自己脑子中，虽然目前看似乎多余，但它
能“武装”你的思想，总有一天会用到。 回想起来，初次接触MATLAB好象是在1996年。当时是因为参加数学建模比赛，其时首先用的还
是Mathematica(还是DOS版本的)，那里MATLAB好象刚出WINDOWNS版本的。后来Mathematica，MATHCAD都用过，
但用了MATLAB后，就不再用另外两个了（汗！）。 要说体会，我可以总结以下几条： 1．多动手写程序、调试 2．善于利用MATLAB的帮助 3．善于向别人学习 4．时间积累 ＊多动手写程序、调试 如
果懒得写程序，调试程序，永远无法提高。我个人认为调试程序更重要。有些朋友可能在一个程序调试几下出不了结果时，就可能喜欢去问别人，我不太赞同这一做
法。其实，凡事往往经过痛苦折磨后，才会让你印象深刻，收益更大。我建议在你觉得用尽你努力后，仍然无法有结果时，才去请教别人。我当初一个程序调试过一
两个星期都有过。在这论坛上，你可以发现不少好的问题，对这些问题，不要光看别人如果解决，也不要光想怎么解决，自己坐下来，动手自己解决一下，那你就会
把不是你的知识变成自己的知识。 ＊善于利用MATLAB的帮助 可以这么说，任何问题都可以在MATLAB的帮助里找到解决的办
法。问题不论大小，都是由更小的问题组成，把大问题化为小问题，小函数，然后再到MATLAB帮助里去找这种小问题，小函数的用法。说实话，MATLAB
里的函数太多，我也经常忘记一些用法，这时HELP就帮忙了。 ＊善于向别人学习 在你解决一个问题后，你可能会发现别人有更简便的方法解决，更强的函数，就时就是你向别人学习的时候。说实话，在这论坛上，我也向bzzz, bainhome等学习借鉴不少。 ＊时间积累 时间长了，积累多了，当然也就有进步了。呵呵，也许再过几年，你会发现原来问题也不是以前想的那么难。而lyrock在这里发的也是“打糊乱说，小儿科”，那时你就已经积累不少了。 我不想说“趁自己年轻，多学点”这种话，想说的是“活到老，学到老”，呵呵，因为我已经不是很年轻了。祝大家共同进步！

redstar · 发表于 2005-8-21 14:52

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 如有疏漏，请指正^_^

redstar · 发表于 2005-8-21 14:53

今天在研学上看到的一个问题 <blockquote>[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 问题：因为不考虑任何摩擦和阻尼，求解该问题的结果应当是一正弦曲线的样子，不会随时间衰减的，而现在衰减了，请教问题出在什么地方？另外图片怎样加上去？我怎么加不上去 </blockquote> 首先说明一个问题 这样写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取得这么长没什么意义，这么简单的方程不如用符号函数 下图蓝色的是减小误差容限，红色的是默认容限
 
<a href="http://www.simwe.com/forum/upload/2005/08/11/83685021.jpg" target="_blank" ><img src="http://www.simwe.com/forum/upload/2005/08/11/83685021.snap.jpg" title="Click to view full a.jpg (1201 X 900)" align="middle" border="0" height="479" width="640"></a> （缩略图，点击图片链接看原图）

redstar · 发表于 2005-8-21 14:53

在正餐开始前,先来点暖场开胃的,有关euler方法的MATLAB实现和比较(忙了一下午,从前自学的东东都忘完了...^_^): function diff_numeric_solve_euler_1 %method1:euler法---------y(n+1)=yn+h*f(xn,yn)--->显式计算公式 xn=0:.1:1; y0=1; %赋初值 y_n_plus1=[]; for i=1:length(xn)-1 y_n_plus1=[y_n_plus1;y0+diff(xn([1:2]))*(y0-2*xn(i)/y0)]; y0=y_n_plus1(end); end y_n_plus1; %method2:后退euler法---------y(n+1)=yn+h*f(x(n+1),y(n+1))--->隐式计算公式 xn1=0:.1:1; y01=1; %赋初值 y_n2=[]; for i=1:length(xn1)-1 y02=y01+diff(xn1([1:2]))*(y01-2*xn1(i)/y01); y_n2=[y_n2;y01+diff(xn1([1:2]))*(y02-2*xn1(i+1)/y02)]; y01=y_n2(end); end %上述两种方法的局部截断误差分别约为h^2*y''(xn)/2,-h^2*y''(xn)/2 %利用中心差分公式的梯形公式可以得到更好的精度 xn2=0:.1:1; y11=1; y_n3=[]; for i=1:length(xn2)-1 y03=y11+diff(xn2([1:2]))*(y11-2*xn2(i)/y11); y_n3=[y_n3;y11+.5*diff(xn2([1:2]))*(y11-2*xn2(i)/y11+y03-2*xn2(i+1)/y03)]; y11=y_n3(end); end %建立预测-校正系统的改进euler公式 xnp=0:.1:1; y21=1; y_n4=[]; for i=1:length(xnp)-1 y04=y21+diff(xnp([1:2]))*(y21-2*xnp(i)/y21); y_n4=[y_n4;y21+.5*diff(xnp([1:2]))*(y21-2*xnp(i)/y21+y04-2*xnp(i+1)/y04)]; y21=y_n4(end); end %精确的解析解 dy=dsolve('Dy=y-2*x/y','y(0)=1','x'); y_exact=subs(dy,'x',{.1:.1:1}); %ode45的解 f=inline('x-2*t/x','t','x'); [t,Y]=ode45(f,[0,1],1); cc=flipdim(Y([end:-4:1]),1)'; y_ode=cc([2:end]); %精度比较 y_compare=[y_n_plus1';y_n2';y_n3';y_n4';y_ode;y_exact]; y_compare_minor=[y_compare(6,:)-y_compare(1,:);... y_compare(6,:)-y_compare(2,:);y_compare(6,:)-y_compare(3,:);... y_compare(6,:)-y_compare(4,:);y_compare(6,:)-y_compare(5,:)]; var_y=var(y_compare_minor'); =============================================================================== 以下是运行结果: y_compare 1.1 1.1918 1.2774 1.3582 1.4351 1.509 1.5803 1.6498 1.7178 1.7848 1.0918 1.1763 1.2546 1.3278 1.3964 1.4609 1.5216 1.5786 1.6321 1.6819 1.0959 1.1841 1.2662 1.3434 1.4164 1.486 1.5525 1.6165 1.6782 1.7379 1.0959 1.1841 1.2662 1.3434 1.4164 1.486 1.5525 1.6165 1.6782 1.7379 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 算出的五种方法的方差为: var_y=0.00025172 0.00023498 3.1078e-006 3.1078e-006 5.1132e-018 显
然可以看出,第一,二种方法的精度比较差,第三,四种有明显改善,且这两种方法实际是一个意思,建立预-校系统的解法算法由于提供了一个预测值,相对来说
有简化.第五种方法,也就是ode45所得的精度最高尺度到了1e-18...简直可怕!一般很少有工程问题是它不能够满足的了,能够看出,解的差别在本
问题里已经很明显了. 最后还有就是有个小问题请教bzzz,怎么我用inline和@做内联函数使ode45调用老是出错,害我要按传统的方法再
写一个外部函数fun_ode,晕倒...你给个建议吧!最好在程序里直接改,我就可以复制粘贴了...the way I go just
described as below:^_^ f=inline('x(2)-2*x(1)/x(2)','t','x'); ================================================================================ 按
照bzzz的建议对ode45的地方做了改动简化,且今天早上把改进的euler方法添加进去了,以上都是单步法的计算公式,euler两步法和线性多步
法的东东就不再介绍,从数值计算的方向考虑道理相同,只是循环次数改为n-2,y值多计算一次.以上方法均满足李普希茨稳定条件,即:y=y(x),可以
存在且唯一.而后退euler法更是无条件稳定的,稳定性和代数精度方面的内容哪位有比较深的研究欢迎发贴给予总结.

账号		自动登录	找回密码
密码			我要加入

[编程技巧] 常微分,偏微分方程的数值解法讨论 zz

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分