2017-04-21 115 views
0

下面是我试过的代码:数值答案 - 微分方程MATLAB

syms y(x) 
Du = diff(y,x); 
ode = diff(y,x,2) - (0.5/(x+1))*diff(y,x)+0.5*x*y == x; 
cond1 = y(1.3) == 0.5; 
cond2 = Du(1.3) == 2; 
conds = [cond1 cond2]; 

uSol(x) = dsolve(ode,conds) 

如果有一个象征性的答案,将工作,但我想没有。任何人都可以告诉我如何得到这个方程的答案?

这就是答案,我得到:

Warning: Explicit solution could not be found. 
> In `dsolve` (line 201) 
    uSol(x) = 
    [ empty sym ] 

+1

对于这类问题,你为什么不使用['钨alpha'](http://www.wolframalpha.com/)? – qbzenker

回答

1

您可以使用ode45或者你可以使用一个数值方法。 0.5 * X * Y + 0.5 -

要将方程分裂成两个一阶方法

U(X)= Y '(x)的

U'(x)= X

无论哪种方式* U /(X + 1)

Y(1.3)= 0.5

U(1.3)= 2

这样做时,自己...

function StackOverflow 
close all 
set(groot,'defaultLineLineWidth',3) 

dx=0.01; 
x = (0:dx:1.3)'; 
N = length(x); 

Y = nan(size(x,1),2); 
Y(N,:) = [0.5, 2]; 

F [email protected](x,y) [y(2) , x-0.5*x.*y(1)+0.5*y(2)./(x+1)]; 

for i=N:-1:2        % calculation loop 
    k_1 = F(x(i), Y(i, :)); 
    k_2 = F(x(i) + 0.5*dx, Y(i,:) + 0.5*dx*k_1); 
    k_3 = F(x(i) + 0.5*dx, Y(i,:) + 0.5*dx*k_2); 
    k_4 = F(x(i) +  dx, Y(i,:) +  dx*k_3); 

    Y(i-1,:) = Y(i,:) + (dx/6)*(k_1+2*k_2+2*k_3+k_4); % main equation 
end 

y = Y(:,1); 
u = Y(:,2); 

figure(1) 
set(gcf, 'units', 'normalized') 

subplot(2,1,1) 
plot(x , y) 
set(gca,'fontsize',14, 'Position',[0.12 0.57 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 
xticklabels({}) 
ylabel('$y(x)$', 'Interpreter', 'latex') 


subplot(2,1,2) 
plot(x , u) 
set(gca,'fontsize',14, 'Position',[0.12 0.125 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 

xlabel('$x$', 'Interpreter', 'latex') 
ylabel('$y''(x)$', 'Interpreter', 'latex') 

end 

,或者让Matlab的为你做...

function StackOverflow2 

%Note: t = xMax - x 
[t,Y] = ode45(@F ,[0 1.3],[0.5; 2]); 

z = flipud([t,Y]); 
x = max(t) - z(:,1); 
y = z(:,2); 
u = z(:,3); 

figure(2) 
subplot(2,1,1) 
plot(x , y) 
set(gca,'fontsize',14, 'Position',[0.12 0.57 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 
xticklabels({}) 
ylabel('$y(x)$', 'Interpreter', 'latex') 


subplot(2,1,2) 
plot(x , u) 
set(gca,'fontsize',14, 'Position',[0.12 0.125 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 

end 

function dydt = F(t,y) 

dydt = [y(2); (1.3-t)-0.5*(1.3-t).*y(1)+0.5*y(2)./((1.3-t)+1)]; 

end 
+0

哦哇..谢谢!我还没有在编程方面进步,所以我不了解它的大部分,代码也没有运行,所以我可能不会找到这个错误。 – Rih

+0

哪一个不运行?什么是错误? – user1543042

0

在最早的Matlab的版本中,函数的定义是不允许里面的脚本。从@user1543042

以答案你应该把这个在一个文件

function StackOverflow2 

%Note: t = xMax - x 
[t,Y] = ode45(@F ,[0 1.3],[0.5; 2]); 

z = flipud([t,Y]); 
x = max(t) - z(:,1); 
y = z(:,2); 
u = z(:,3); 

figure(2) 
subplot(2,1,1) 
plot(x , y) 
set(gca,'fontsize',14, 'Position',[0.12 0.57 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 
xticklabels({}) 
ylabel('$y(x)$', 'Interpreter', 'latex') 


subplot(2,1,2) 
plot(x , u) 
set(gca,'fontsize',14, 'Position',[0.12 0.125 0.85 0.40]) 
xlim([min(x), max(x)]) 
xticks(linspace(min(x), max(x), 5)) 

end 

而这在另一个文件:

function dydt = F(t,y) 

dydt = [y(2); (1.3-t)-0.5*(1.3-t).*y(1)+0.5*y(2)./((1.3-t)+1)]; 

end