求解微分方程后为何画不出求解结果

semeng

function chapter7()
clc
clear all
close all
global tc w
tc=0.1;w=6.28;
final_time=1;
a=[0,final_time];
y0=0;
[t,y]= ode45(@eq4,a,y0);
plot(t,y) % 本想 plot (t,y,t,F,':')
%% plot (t,y,t,F,':') % 这个地方就出错了，为什么找不到求解方程的F值那？可是
%% 运行时去掉F表达式后面的分号，在command窗口有F的求解结果，怎么就画不出来哪？
%% 希望明白人给个明白的解释，谢谢
xlabel('t'),ylabel('y(t)')
axis tight
end

function ydot=eq4(t,y)
        global tc w
        F=t*exp(-t/tc)*cos(w*t);
        ydot=5*(F-y);
end

[ 本帖最后由 eight 于 2008-5-7 11:31 编辑 ]

ch_j1985

原帖由 semeng 于 2008-5-7 09:49 发表

                               
登录/注册后可看大图

function chapter7()
clc
clear all
close all
global tc w
tc=0.1;w=6.28;
%% plot (t,y,t,F,': ...

虽然command窗口有F的求解结果，但是Workspace里面没有F，所以才会报错！
之所以会这样，我想应该是函数体可以调用主程序的全局变量，但是主程序不可以调用函数体的中间变量的缘故吧！

代码应改为这样：

function chapter7()
clc
clear all
close all
global tc w
tc=0.1;w=6.28;
final_time=1;
a=[0,final_time];
y0=0;
[t,y]= ode45(@eq4,a,y0);
F=t.*exp(-t/tc).*cos(w*t);
plot (t,y,t,F,':')
xlabel('t'),ylabel('y(t)')
axis tight
end

[ 本帖最后由 ch_j1985 于 2008-5-7 10:29 编辑 ]

semeng

首先谢谢你的解答，我现在已经改过来了。但是还有一个问题需要您继续回答
请问，F的表达式实际上在方程（eq4）里面已经出现过了
而且方程也能求出F的每一个数值
为什么就不能提取
如果F表达式很复杂，重新编排，整理（我是说变成矩阵运算）岂不是相当麻烦，
还有没有不用再输入一次表达式的方法呢？
谢谢您的关注

[ 本帖最后由 eight 于 2008-5-7 10:55 编辑 ]

wwbeyondww

把F也定义成全局变量行不行啊?

semeng

no
It has some errors
please, try it and explain it

wwbeyondww

这样写好像可以,呵呵

function chapter7()
clc
clear all
close all
global tc w
tc=0.1;w=6.28;
final_time=1;
a=[0,final_time];
y0=0;
[t,y]= ode45(@eq4,a,y0);
[F,a]=eq4(t,y);
% plot(t,y) %

plot (t,y,t,F,':')
%% plot (t,y,t,F,':') % 这个地方就出错了，为什么找不到求解方程的F值那？可是
%% 运行时去掉F表达式后面的分号，在command窗口有F的求解结果，怎么就画不出来哪？
%% 希望明白人给个明白的解释，谢谢
xlabel('t'),ylabel('y(t)')
axis tight
end

function [F,ydot]=eq4(t,y)
        global tc w
        F=t.*exp(-t/tc).*cos(w*t);
        ydot=5*(F-y);
end

eight

原帖由 wwbeyondww 于 2008-5-7 11:00 发表

                               
登录/注册后可看大图

这样写好像可以,呵呵

function chapter7()
clc
clear all
close all
global tc w
tc=0.1;w=6.28;
final_time=1;
a=[0,final_time];
y0=0;
[t,y]= ode45(@eq4,a,y0);
[F,a]=eq4(t,y);
% plot(t,y) %
...

有出错提示吗？没有的话，自己调试一下吧

semeng

很多错误提示
我都不知道怎么改了，
学识不深，请直言

ch_j1985

原帖由 semeng 于 2008-5-7 15:30 发表

                               
登录/注册后可看大图

很多错误提示
我都不知道怎么改了，
学识不深，请直言

按照5楼的代码，编程画出的图如下（不知是不是你想要的）：

wwbeyondww

又没改啥,只是多加了个函数返回的参数而已啊