发信人: hahn (有奇☆闭关造文中), 信区: Matlab
标 题: [合集] Matlab解微分方程问题求助
发信站: 哈工大紫丁香 (Tue Nov 14 07:34:10 2006), 站内
────────────────────────────────────────
COWOO (蔚蓝的风) 于 (Mon Jun 5 21:17:36 2006) 说道:
高中同学让帮忙算个方程画个图，好就没摸Matlab了，微分方程也不太熟，不知道什么问
题。
同学给的题目如下：
边界条件:
r=0,ds/dr=0
r=R,s=5
第二个边界条件也可以是ds/dr＝无穷大 其中m,n,h都是常数 你可以取的 你先随便取 看
看画出来的图怎么样.
微分方程：
<img src="http://blog.cs.hit.edu.cn/function.png"/>
其中r的范围是0到0.00225，s的范围是0到5
我的解法：
用龙格库塔方法，把这个方程化成方程组，并在写成一个叫srfun.m的函数
function zhangfan=srfun(r,x);
global m
global n
global h
zhangfan=[x(2);m*x(1)*exp(-n*r)+h*x(1)-2*x(2)/r]
然后在matlab中调用：
global m;
m=1;
global n;
n=1;
global h;
h=1;
[r,s]=ode23('srfun',[0 0.00225],[2 0]);
结果算出来s除了第一组值是［2 0］，其余全部是NaN。
我认为可能出错的地方：
1. 微分方程化错了，好久没动微分方程了，化错了大家不要笑话我，呵呵
2. 初值［2 0］选的不对。请高手指点这个问题怎么选这个参数？
3. 还有一个可能，就是r的范围太小，数量级相差太大，但是我把［0 0.00225］换成
［0 2.25］ 还是一样效果，我不知所措了。。。
望高手指点一二。
────────────────────────────────────────
zjliu (秋天的萝卜) 于 (Mon Jun 5 21:53:49 2006) 说道:
Example：
% 二分法+打靶法解微分方程
% 方程:
% diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
% where diff is a difference for r.
% 边界条件:
%
% r=0,ds/dr=0
% r=R,s=5
clc;clear;close all;
m=1;
n=1;
h=1;
R=5;
fun=inline('[s(2);m*s(1)*exp(-n*r)+h*s(1)-2*s(2)/r]',...
'r','s','flag','m','n','h');
s_end=5;
sp=1;
rsa=0.1;
rsb=1;
d=1;
[r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
[r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
figure;
subplot(121);
plot(r1,s1(:,1));hold on;
plot(r2,s2(:,1),'r');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),', rsb=',num2str(rsb)]);
subplot(122);
pa=plot(r1,s1(:,1));hold on;
pb=plot(r2,s2(:,1),'r');
pc=plot(r2,s2(:,1)*0.6,'k');
legend('for rsa','for rsb',0);
xlim([0,R]);
tg=title(['rsa=',num2str(rsa),', rsb=',num2str(rsb)]);
while abs(s1(end,1)-s2(end,1))>1e-4; % 控制精度
[r1,s1]=ode45(fun,[eps,R],[rsa;0],[],m,n,h);
[r2,s2]=ode45(fun,[eps,R],[rsb;0],[],m,n,h);
rs=[rsa+rsb]/2;
[r,s]=ode45(fun,[eps,R],[rs;0],[],m,n,h);
if s(end,1)>5; % 二分法部分
rsb=rs;
else
rsa=rs;
end
set(pc,'XData',r,'YData',s(:,1));
set(tg,'string',['rs=',num2str(rs)]);
pause(0.2);
end
title(['This program is over!'])
【 在 COWOO (蔚蓝的风) 的大作中提到: 】
: 高中同学让帮忙算个方程画个图，好就没摸Matlab了，微分方程也不太熟，不知道什么问
: 题。
: 同学给的题目如下：
: ...................
────────────────────────────────────────
COWOO (蔚蓝的风) 于 (Mon Jun 5 22:18:50 2006) 说道:
这么快！真是太感谢了。
赞一个！
【 在 zjliu (秋天的萝卜) 的大作中提到: 】
: Example：
: % 二分法+打靶法解微分方程
: % 方程:
: % diff(s,2)+2*diff(s,1)=m*s*exp(-n*r)+h*s,
: % where diff is a difference for r.
: % 边界条件:
: %
: % r=0,ds/dr=0
: ...................
────────────────────────────────────────
feifeifool (爱芳) 于 (Tue Jun 6 09:52:48 2006) 说道:
去处m*s*exp(-r)*s,是个欧拉方程,是有解析解的.
不知道这个变系数的项,能不能用常数变异法求解出来,有理论学过常微分的同学能不能执教一下.
【 在 COWOO (蔚蓝的风) 的大作中提到: 】
: 高中同学让帮忙算个方程画个图，好就没摸Matlab了，微分方程也不太熟，不知道什么问
: 题。
: 同学给的题目如下：
: ...................
────────────────────────────────────────
zjliu (秋天的萝卜) 于 (Tue Jun 6 10:13:19 2006) 说道:
你说的欧拉方程，计算它的本征值方程就可以计算了
【 在 feifeifool (爱芳) 的大作中提到: 】
: 去处m*s*exp(-r)*s,是个欧拉方程,是有解析解的.
: 不知道这个变系数的项,能不能用常数变异法求解出来,有理论学过常微分的同学能不能执教一下.