solve('x^4+x+c=0','x')

clear
clc

%输入平台系统参数
display('输入平台系统参数');
% d=input('太阳光吸收涂层宽度d(m):');
% D=input('玻璃管的直径D(m):');
% D_tro=input('反射镜宽度(m):');
% L=input('平台系统长度(m):');
% E_p=input('太阳光选择吸收涂层的辐射率:');
% E_g=input('玻璃的辐射率:');
% E_tro=input('反射镜的辐射率:');
% E_rear=input('冷却管背部的辐射率:');
% a_rear=input('冷却管背部的吸收率:');
d=0.06;
D=0.12;
D_tro=3;
L=6;
E_p=0.05;
E_g=0.2;
E_tro=0.2;
E_rear=0.2;
a_rear=0.8;

%输入周围环境参数
display('输入周围环境参数')
% v=input('周围环境风速(m/s):');
% Ta=input('环境空气温度Ta(K)：');
% Tc=input('冷却管温度Tc(K):');
v=1.5;
Ta=300;
Tc=323;

%常数
%斯蒂芬－波尔兹曼常数
sigma=5.667e-8;  % W/(m2.K4)



%定义玻璃管温度Tg为变量
syms Tg Tp


%%%太阳能热电系统热损耗系数如下%%%
%（1）太阳光选择吸收涂层与玻璃管下部之间辐射：
hr_sg=sigma*(Tp^2+Tg^2)*(Tp+Tg)/(1/E_p+1/E_g-1);
Ap=d*L;
R1=1/(hr_sg*Ap);
%（2）玻璃管上部与冷却管之间的辐射：
hr_rg=sigma*(Tc^2+Tg^2)*(Tc+Tg)/(E_rear+E_g-1);
R2=1/(hr_rg*Ap);
%（3）玻璃管下部与反射镜之间辐射：
hr_gtro=sigma*(Tg^2+Ta^2)*(Tg+Ta)/(E_g+E_tro-1);
Ag=pi*D*L;
R3=1/(hr_gtro*0.5*Ag);
%（4）玻璃管下部与反射镜之间空气对流：
hw=5.7+3.8*v;  % W/(m2.K)
R4=1/(hw*0.5*Ag);
%(5)反射镜与周围环境空气对流：
%略
%（6）玻璃管上部对天空的辐射：
hr_gsky=sigma*(Tg^2+Ta^2)*(Tg+Ta)*E_g;
R6=1/(hr_gsky*0.5*Ag);
%（7）反射镜对天空的辐射：
%略
%（8）玻璃管上部与周围环境空气对流：
R8=1/(hw*0.5*Ag);

%求解Tg
f=(Tp-Tg)/R1+(Tc-Tg)/R2-(Tg-Ta)*(1/R3+1/R4+1/R6+1/R8);

%求数值解
tp = 450:650; %因为上面f里面的Tp还有用，所以这里不能同名
tg = zeros(1, length(tp));

for i = 1:length(tp)
    tf = subs(f, Tp, tp(i)); %给Tp代入值
    fHandle = inline(char(tf), 'Tg'); %将tf视为函数句柄，'Tg'为自变量
   
    %在Tg=400附近的实数域找解，如果找不到会返回NaN
    tg(i) = fzero(fHandle, 400);
   
    if tg(i)<=400
        tg(i) = 0; %如果没找到大于400的解，那么就置为0表示没找到
    end
end

plot(tp, tg);

for i=1:length(y)
    if isreal(y(i)) && y(i)>400 %判断y(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f\n', y(i));
    end
end

syms x y; %注意用到的符号都要声明一下
a=y+1;
b=y-1;

f=a*b+a+b*x+y+6*x+2;

%解方程f=0
equation = [char(f) ' = 0'];
A=solve(equation, 'y');

%上面解得的结果是
%A =
%  - x/2 - (x^2 - 16*x - 4)^(1/2)/2 - 1
%   (x^2 - 16*x - 4)^(1/2)/2 - x/2 - 1

%对于符号式子，使用subs函数代入数值

%注意上面的方程因为有开方，所以不一定能解出实数解，这里直接取一段有实数解的
x = 20:80;
y1 = zeros(1, length(x)); %习惯上循环前先给定存储空间
y2 = zeros(1, length(x));
for i=1:length(x)
    y1(i) = subs(A(1), x(i));
    y2(i) = subs(A(2), x(i));
end

plot(x, y1, 'b', x, y2, 'r');

x = 20:80;
Y= subs(A, x);
plot(x, Y(1, :), 'b', x, Y(2, :), 'r');

2楼: Originally posted by sudo at 2011-10-12 10:12:19:
你不用每次都放那么多分上来的...

(1)使用isreal函数和循环
[code]
for i=1:length(y)
    if isreal(y(i)) && y(i)>400 %判断y(i)是否是实数，并且是否大于400
        fprintf(1, '%.9f ...