MATLAB微分方程的求解,请解释一下这个程序: function [c,f,s]=Pdepe(x,t,u,DuDx) %PDE方程函数 global a

MATLAB微分方程的求解,请解释一下这个程序: function [c,f,s]=Pdepe(x,t,u,DuDx) %PDE方程函数 global a,第1张

pdefun为使用者提供的pde函数文件,函数格式为[c,f,s]=pdefun(x,t,u,dudx)使用者仅需提供偏微分方程中的系数向量。c,f和s均为行向量,而向量c即为矩阵c的对角线元素。

显式前向欧拉法

源程序:

function [u,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N)

%解方程 A u_xx = u_t , 0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M : x 轴的等分段数

% N : t 轴的等分段数

dx = xf/M; x = [0:M]dx;

dt = T/N; t = [0:N]'dt;

for i= 1:M + 1

u(i,1) =it0(x(i));

end

for j = 1:N + 1

u([1 M + 1],j)= [bx0(t(j)); bxf(t(j))];

end

r = Adt/dx/dx, r1 = 1 - 2r;

if(r>05)

disp('r>05,unstability');

end

for j = 1:N

for i = 2:M

u(i,j+1) =r(u(i + 1,j) + u(i-1,j)) + r1u(i,j); %(923)

end

end

u=u';

%------------------------------------

在MATLAB中编写脚本文件:

A = 05; %方程系数

it0 = inline('sin(pix)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 80; T = 01; N = 100;

[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);

figure(1),clf,mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

title('r>05')

M=50;

[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);

figure(2),clf,mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

title('r<05')

%-----------------

隐式后向欧拉法

源程序:

function [u,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N)

%解方程 A1 u_xx = u_t , 0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M : x 轴的等分段数

% N : t 轴的等分段数

dx = xf/M; x = [0:M]dx;

dt = T/N; t = [0:N]'dt;

for i = 1:M + 1

u(i,1) =it0(x(i));

end

for j = 1:N + 1

u([1 M + 1],j)= [bx0(t(j)); bxf(t(j))];

end

r = Adt/dx/dx;

r2 = 1 + 2r;

for i = 1:M - 1

P(i,i) = r2; %构造929的矩阵

if i > 1

P(i - 1,i)= -r; P(i,i - 1) = -r;

end

end

for j = 2:N + 1

b = [ru(1,j);zeros(M - 3,1); ru(M + 1,j)] + u(2:M,j - 1); %Eq(929)

u(2:M,j) =linsolve(P,b);

end

u=u';

%---------------------

在MATLAB中编写脚本文件:

A= 05; %方程系数

it0 = inline('sin(pix)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 50; T = 01; N = 100;

[u1,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N);

mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

%--------------

Grank-Nicholson方法

源程序:

function [u,x,t] =Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N)

%解方程 A u_xx = u_t , 0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M : x 轴的等分段数

% N : t 轴的等分段数

dx = xf/M; x = [0:M]dx;

dt = T/N; t = [0:N]'dt;

for i = 1:M + 1

u(i,1) =it0(x(i));

end

for n = 1:N + 1

u([1 M + 1],n)= [bx0(t(n)); bxf(t(n))];

end

r = Adt/dx/dx;

r1 = 2(1 + r); r2 = 2(1 - r);

for i = 1:M - 1

P(i,i) = r1; %(9217)

Q(i,i)=r2;

if i > 1

P(i - 1,i)= -r; P(i,i - 1) = -r; %(9217)等式左边矩阵

Q(i - 1,i)= r; Q(i,i - 1) = r;%(9217)等式右边矩阵

end

end

for k = 2:N + 1

b=Qu(2:M,k-1)+[r(u(1,k)+u(1,k-1));zeros(M-2,1)];

u(2:M,k) =linsolve(P,b); %(9217)

end

u=u';

%-----------

在MATLAB中编写脚本文件:

A= 05; %方程系数

it0 = inline('sin(pix)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 25; T = 01; N = 100;

[u1,x,t] = Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N);

mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

请问用matlab怎么解

dx/dt=2y^2-2

dy/dt=2yz-x-1

dz/zt=2z^2-2y-x0=[1

0

2];

%微分变量区间

tspan=[3

0];

%注意这里必须写成[3

0]

以上就是关于MATLAB微分方程的求解,请解释一下这个程序: function [c,f,s]=Pdepe(x,t,u,DuDx) %PDE方程函数 global a全部的内容,包括:MATLAB微分方程的求解,请解释一下这个程序: function [c,f,s]=Pdepe(x,t,u,DuDx) %PDE方程函数 global a、急求!!! 大学数学,用matlab解决问题,题目是一维抛物型偏微分方程差分解法、求 MATLAB程序 用有限元方法解二阶偏微分方程,初值问题等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址:https://54852.com/zz/9995562.html

(0)
打赏 微信扫一扫微信扫一扫 支付宝扫一扫支付宝扫一扫
上一篇 2023-05-04
下一篇2023-05-04

发表评论

登录后才能评论

评论列表(0条)

    保存