如何用matlab求解微分方程组的数值解

Posted 2023-04-22

在Matlab下输入：edit
zhidao_feiying.m，然后将下面两行百分号之间的内容，复制进去，保存

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function
y=zhidao_feiying(t,x)

a=1;

b=2;

c=3;

d=4;

%%比方说

%f=a*x+y-z;

%g=b*sin(x*y)-c*cos(z);

%h=d*y-a*x;

%%注意x用x(1)代，y用x(2)代,z用x(3)代

f=a*x(1)+x(2)-x(3);

g=b*sin(x(1)*x(2))-c*cos(x(3));

h=d*x(2)-a*x(1);

y=[f;g;h];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

在Matlab下面输入：

t_end=5;

x0=[1;1;1];

[t,x]=ode45('zhidao_feiying',[0,t_end],x0);

plot(t,x)

legend('x','y','z')

上面只是固定d的情况，

你如果想做出随d的变化，估计比较麻烦一些，一方面ode45并不是等不长的。

另一方面，参数不好弄。 参考技术A 在matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
y=zhidao_rk4_tangwei0206302(t,x)
x1=x(1);
x2=x(2);
x3=x(3);
y=[-x2*(x3-1);
x1*(x3-3);
x1*x2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
在matlab下面输入：
t_end=200;
x0=[1;sqrt(6);0];
[t,x]=ode45('zhidao_rk4_tangwei0206302',[0,t_end],x0);
plot3(x(:,1),x(:,2),x(:,3));
xlabel('x1');
ylabel('x2');
zlabel('x3');

微分方程数值解用matlab练习

实习题1：

将xx化成一阶方程组，并用Euler法和改进的Euler法求解，步长h = 0.1, 0.05。

画出Euler法，Euler改进法以及精确解的数值结果图形

clear all;
clc;
h1 = 0.05;
len = 1 + 1 / h1;
u = zeros(1, len);
v = zeros(1, len);
t = zeros(1, len);
exact = zeros(1, len);

u(1) = 0;
v(1) = 1;
t(1) = 0;
exact(1) = 0;
for i = 1:len-1
    t(i + 1) = i * h1;
    u(i + 1) = u(i) + h1 * v(i);
    v(i + 1) = v(i) - h1 * u(i);
    exact(i + 1) = sin(h1 * i);
end

u1 = zeros(1, len);
v1 = zeros(1, len);
u1(1) = 0;
v1(1) = 1;
for i = 1:len-1
    u1(i + 1) = ((1 + h1 * h1 / 4) * u(i) + h1 * v(i)) / (1 - h1 * h1 / 4);
    v1(i + 1) = v1(i) + (u1(i) + u1(i + 1)) * h1 / 2;
end
    
plot(t, u, ‘r*‘, ‘markersize‘, 10);
hold on
plot(t, u1, ‘r.‘, ‘markersize‘, 20);
hold on
plot(t, exact);
grid on

View Code