你能寄给我的照片第二学位方程组和边界条件?
PDE数字解决方案的曲面图
1视图(30天)
显示旧的评论
你好,我是解决2 PDE的系统,每个函数取决于3个变量(半径、角度和时间)。我使用显式差分格式。由于我有两个3 d矩阵(一个用于每个函数)。想象的结果,我想情节表面情节与固定每个函数t(时间)。例如:在一个时刻t = 0.5 u的表面(:,:0.5):轴半径、轴角和z轴将在这个点的函数值(x, y)。会很高兴的任何建议。谢谢。
这是代码,如果是很难理解我的英语。
如果真正的
%的代码
clc;
清晰的所有;
为r %网格
n = 100;
r_min = 0.01;
r_max = 1;
hr = (r_max-r_min) / (n - 1)
r = r_min:人力资源:r_max
nr = max(大小(r))
%网格为θ
m = 100;
th_min = 0;
th_max = 2 *π;
hth = (th_max-th_min) / (m - 1)
θ= th_min: hth: th_max
n = max(大小(θ))
% t网格
l = 10000;
t_min = 0;
t_max = 1;
ht = (t_max-t_min) / (l - 1)
时间= t_min: ht: t_max
元= max(大小(时间))
u = 0 (nr、n、nt);
v = 0 (nr、n、nt);
% Inititalization
情况= 0;
为我= 1:nr
为j = 1: n
如果(r(我)< r0)
情况= 0;
elseif((r (i) > = r0) & & (r (i) < r1))
情况= (r(我)r0) / (r1-r0);
elseif((r (i) > = r1) & & (r (i) < r2))
情况= 1;
elseif((r (i) > = r2) & & (r (i) < r3))
情况= (r(我)r3) / (r2-r3);
其他的
情况= 0;
结束
1 u (i, j) =情况;
结束
结束
u (:: 1)
v0 = 0;
为我= 1:nr
为j = 1: n
如果(r(我)< r1)
= 1;
elseif((r (i) > = r1) & & (r (i) < r3))
v0 = (r(我)r3) / (r1-r3);
其他的
v0 = 0;
结束
(i, j, 1) = v0;
结束
结束
v (:: 1)
为t = 1: nt-1
为我= 2:nr-1
为j = 2: nth-1
%的衍生品
d2udr2 = (u (i + 1, j, t) 2 * u (i, j, t) +(张,j, t)) /小时^ 2;
dudr = (u (i + 1, j, t) - u(张,j, t)) /(2 *人力资源);
d2udth2 = (u (i, j + 1, t) 2 * u (i, j, t) + (i, j - 1 t)) / hth ^ 2;
d2vdr2 = (v (i + 1, j, t) 2 * v (i, j, t) + v(张,j, t)) /小时^ 2;
dvdr = (v (i + 1, j, t) - v (i, j, t) /人力资源;
d2vdth2 = (v (i, j + 1, t) 2 * v (i, j, t) + v (i, j - 1 t)) / hth ^ 2;
%中心节点
u (i, j, t + 1) = u (i, j, t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u (i, j, t) * (1 u (i, j, t) - c * v (i, j, t)));
v (i, j, t + 1) = v (i, j, t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v (i, j, t) * (1 - c * u (i, j, t) - b * v (i, j, t)));
结束
%为左边界衍生品
d2udr2 = (u (i + 1, 1 t) 2 * u(我1 t) + u(1张,t)) /小时^ 2;
dudr = (u (i + 1, 1 t) - u(1张,t) /(2 *人力资源);
d2udth2 = (u (2, t) 2 * u(我1 t) + u (n, t)) / hth ^ 2;
d2vdr2 = (v (i + 1, 1 t) 2 * v(我1 t) + v(张1 t) /小时^ 2;
dvdr = (v (i + 1, 1 t) - v(张1 t) /(2 *人力资源);
d2vdth2 = (v(我2 t) 2 * v(我1 t) + v (n, t) / hth ^ 2;
uleft = u(我1 t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u(我1 t) * (1 u(我1 t) - c * v(我1 t)));
vleft = v(我1 t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v(我1 t) * (1 - c * u(我1 t) - b * v(我1 t)));
%衍生品右边界
d2udr2 = (u (i + 1, n, t) 2 * u (n, t) + u(张、n t)) /小时^ 2;
dudr = (u (i + 1, n, t) - u(张,n, t)) /(2 *人力资源);
d2udth2 = (u(我1 t) 2 * u (n, t) + u(我nth-1 t)) / hth ^ 2;
d2vdr2 = (v (i + 1, n, t) 2 * v (n, t) + v(张,n, t)) /小时^ 2;
dvdr = (v (i + 1, n, t) - v(张,n, t)) /(2 *人力资源);
d2vdth2 = (v(我1 t) 2 * v (n, t) + v(我nth-1 t) / hth ^ 2;
uright = u (n, t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u (n, t) * (1 u (n, t) - c * v (n, t)));
vright = v (n, t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v (n, t) * (1 - c * u (n, t) - b * v (n, t)));
%填补边界为θ
u(我1 t + 1) = uleft;
v(我1 t + 1) = vleft;
u(我n t + 1) = uright;
v(我n t + 1) = vright;
结束
为r %边界
为j = 1: n
u (1 j t + 1) = u (2 j t + 1);
u (nr, j, t + 1) = u (nr-1 j t + 1);
v (t + 1 1 j) = v (t + 1 2 j);
v (nr, j, t + 1) = v (t + 1 nr-1 j);
结束
结束
结束
冲浪(r,θ,?)
接受的答案
优素福Khmou
2013年12月14日
编辑:优素福Khmou
2013年12月14日
Danila,我试着给你检查解决方案,第一件事就是r0, r1, r2和r3不予提供,加上一个变量μ没有定义,我相信μ是一个内部变量名,所以我把它改为μ的区别,这是改变了版本的解决方案。关于最后一个问题你可以简单地使用
冲浪(r,θ,u (:,:, n))
第一即时n = 1, n = 2, .....n =结束:
clc;明确所有;
n = 100;
r_min = 0.01;
r0 = r_min;
r_max = 1;
r3 = r_max;
r1 = r0 + 0.2;
r2 = r1 + 0.2;
hr = (r_max-r_min) / (n - 1);
r = r_min:人力资源:r_max;
nr = max(大小(r));
%网格为θ
m = 100;
th_min = 0;
th_max = 2 *π;
hth = (th_max-th_min) / (m - 1);
θ= th_min: hth: th_max;
n = max(大小(θ));
% t网格
l = 100;
t_min = 0;
t_max = 1;
ht = (t_max-t_min) / (l - 1);
时间= t_min: ht: t_max;
元= max(大小(时间));
u = 0 (nr、n、nt);
v = 0 (nr、n、nt);
% Inititalization
情况= 0;
为我= 1:nr
为j = 1: n
如果r (r(我)< r0 | | (i) > r3)
情况= 0;
elseif((r (i) > = r0) & & (r (i) < r1))
情况= (r(我)r0) / (r1-r0);
elseif((r (i) > = r1) & & (r (i) < r2))
情况= 1;
elseif((r (i) > = r2) & & (r (i) < r3))
情况= (r(我)r3) / (r2-r3);
结束
1 u (i, j) =情况;
结束
结束
u (:: 1);
v0 = 0;
为我= 1:nr
为j = 1: n
如果(r(我)< r1)
= 1;
elseif((r (i) > = r1) & & (r (i) < r3))
v0 = (r(我)r3) / (r1-r3);
其他的
v0 = 0;
结束
(i, j, 1) = v0;
结束
结束
v (:: 1);
μ= 1;
为t = 1: nt-1
为我= 2:nr-1
为j = 2: nth-1
%的衍生品
d2udr2 = (u (i + 1, j, t) 2 * u (i, j, t) +(张,j, t)) / (hr ^ 2);
dudr = (u (i + 1, j, t) - u(张,j, t)) /(2 *人力资源);
d2udth2 = (u (i, j + 1, t) 2 * u (i, j, t) + (i, j - 1 t)) / (hth ^ 2);
d2vdr2 = (v (i + 1, j, t) 2 * v (i, j, t) + v(张,j, t)) / (hr ^ 2);
dvdr = (v (i + 1, j, t) - v (i, j, t) /人力资源;
d2vdth2 = (v (i, j + 1, t) 2 * v (i, j, t) + v (i, j - 1 t)) / (hth ^ 2);
%中心节点
u (i, j, t + 1) = u (i, j, t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u (i, j, t) * (1 u (i, j, t) - c * v (i, j, t)));
v (i, j, t + 1) = v (i, j, t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v (i, j, t) * (1 - c * u (i, j, t) - b * v (i, j, t)));
结束
%为左边界衍生品
d2udr2 = (u (i + 1, 1 t) 2 * u(我1 t) + u(1张,t)) /小时^ 2;
dudr = (u (i + 1, 1 t) - u(1张,t) /(2 *人力资源);
d2udth2 = (u (2, t) 2 * u(我1 t) + u (n, t)) / hth ^ 2;
d2vdr2 = (v (i + 1, 1 t) 2 * v(我1 t) + v(张1 t) /小时^ 2;
dvdr = (v (i + 1, 1 t) - v(张1 t) /(2 *人力资源);
d2vdth2 = (v(我2 t) 2 * v(我1 t) + v (n, t) / hth ^ 2;
uleft = u(我1 t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u(我1 t) * (1 u(我1 t) - c * v(我1 t)));
vleft = v(我1 t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v(我1 t) * (1 - c * u(我1 t) - b * v(我1 t)));
%衍生品右边界
d2udr2 = (u (i + 1, n, t) 2 * u (n, t) + u(张、n t)) /小时^ 2;
dudr = (u (i + 1, n, t) - u(张,n, t)) /(2 *人力资源);
d2udth2 = (u(我1 t) 2 * u (n, t) + u(我nth-1 t)) / hth ^ 2;
d2vdr2 = (v (i + 1, n, t) 2 * v (n, t) + v(张,n, t)) /小时^ 2;
dvdr = (v (i + 1, n, t) - v(张,n, t)) /(2 *人力资源);
d2vdth2 = (v(我1 t) 2 * v (n, t) + v(我nth-1 t) / hth ^ 2;
uright = u (n, t) + ht * (d2udr2 + (1 / r (i)) * dudr + (1 / r(我)^ 2)* d2udth2);% + u (n, t) * (1 u (n, t) - c * v (n, t)));
vright = v (n, t) + ht *μ* (d2vdr2 + (1 / r (i)) * dvdr + (1 / r(我)^ 2)* d2vdth2);% + * v (n, t) * (1 - c * u (n, t) - b * v (n, t)));
%填补边界为θ
u(我1 t + 1) = uleft;
v(我1 t + 1) = vleft;
u(我n t + 1) = uright;
v(我n t + 1) = vright;
结束
为r %边界
为j = 1: n
u (1 j t + 1) = u (2 j t + 1);
u (nr, j, t + 1) = u (nr-1 j t + 1);
v (t + 1 1 j) = v (t + 1 2 j);
v (nr, j, t + 1) = v (t + 1 nr-1 j);
结束
结束
图,冲浪(r,θ,u(:,:, 3)),阴影插值函数、标题(“U_ {3}”);
图,冲浪(r,θ,v(:,:, 3)),阴影插值函数、标题(“V_ {3}”);
