matlab 微分方程组求解
这个方程满足李普希兹条件，因此，解存在唯一并且可以唯一延拓到边界，应用lax等价定理可以知道，向前欧拉法具有二阶的收敛速度……所以可以尝试用向前欧拉法编写：
我刚学matlab~写的程序一般，希望lz见谅：
这个程序，需要lz选择步长h（我自己试验之后，取h=0.0001的效果比较好），和x_0=1时，y_0的初值（最好取绝对值大于等于2的数）~
% Euler method: numerical method for different equation
% Initialize the initial value of y_0
y_0=input('Initialize y_0=');
% Choose some bandwith h
h=input('Bandwith h=');
% x from x_0=1 to 5: x_up
x_up=1:h:5; % initialize x_up
num=length(x_up); % the number of x that we calculate
y_up=zeros(size(x_up));
y_up(1)=y_0;
ii=1;
while ii<num
y_up(ii+1)=y_up(ii)+((0.3*(1-x_up(ii)^2)*y_up(ii)-x_up(ii))/y_up(ii))*h;
ii=ii+1;
end
% x from -3 to x_0=1: x_down
x_down=-3:h:1;
num=length(x_down);
y_down=zeros(size(x_down));
y_down(num)=y_0;
ii=num;
while ii>1
y_down(ii-1)=y_down(ii)-((0.3*(1-x_down(ii)^2)*y_down(ii)-x_down(ii))/y_down(ii))*h ;
ii=ii-1;
end
x=[x_down,x_up];
y=[y_down,y_up];
while ii<=num
if y(ii)>4 || y(ii)<-3
x(ii)=1;
y(ii)=3;
end
ii=ii+1;
end
plot(x,y,'b+');%可以更加需要修改plot选项
非常感谢您的回答。真是不好意思，深夜你还为我的问题而操心。
我想要的效果如图所示，应该是一个极限环。。但是这需要对初始y值不断的寻找，还请你能帮我解答一下。。
再次感谢！！！
重新写的代码，以微分方程组为基础，用欧拉中心差分近似的结果（x_i+1-x_i-1/2h 代替x在t_i的导数，有二阶的误差，因此，收敛速度是o（h），之前那个有点问题，其实是o（1）的速度）：
% Euler centric different method
% initialize y_0
y_0=input('y_0=');
% select some bandwidth h
h=input('h=');
% select the final point tf
tf=input('tf=');
t=0:h:tf;
num=length(t);
x=zeros(size(t));
y=zeros(size(t));
x(1)=1;
y(1)=y_0;
x(2)=x(1)+y(1)*h;
y(2)=y(1)+(0.3*(1-x(1)^2)*y(1)-x(1))*h;
for ii=3:num
x(ii)=x(ii-2)+y(ii-1)*2*h;
y(ii)=y(ii-2)+(0.3*(1-x(ii-1)^2)*y(ii-1)-x(ii-1))*h*2;
end
plot(x,y);
axis([-3 5 -3 4])
---------------------------------------以上放在一个新建的m-file里面(File\new\script)，名字是eulercdf.m
然后在命令窗口中运行：
>>eulercdf
y_0=0.1
h=0.001
tf=30
>>x1=x;
>>y1=y;
>>eulercdf
y_0=2.5
h=0.001
tf=10
>>x2=x;
>>y2=y;
>>plot(x1,y1,x2,y2);
-----------------------------------------------------------然后得到下面的画图结果：
我先看看有没有大神回答，如果没有我再献丑吧
大神，问题如上追问，请赐教。。。
还好我没献丑……楼上太强了
