现代科学运算-matlab语言与应用
东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解(第四版)》
代码可在matlab r2016b 运行。
02 matlab语言程序设计基础
02.01 数据结构
变量,字母开头,区分大小写。
保留常量
eps, i, j, pi, NaN, Inf, i = sqrt(-1)
lastwarn, lasterr
数值型,双精度 64位,8字节
single(), 32位单精度。
符号型, sym(A),常用于公式推导和求解析解。
变量声明 syms var_list var_props
显示符号变量的任何精度:vpa(A), vpa(A, n)
默认精度–32位10进制有效数字
显示符号变量的属性 assumptions()
设置符号变量的属性 assume(), assumeAlso()
syms x real; assume(x >= -1); assumeAlso(x < 5);
例2-1 符号型数值与双精度数值的区别
1/3 的双精度存储
双精度存储 – 0.3333…33, 15位有限数字
符号型存储 – 1/3:sym(1/3),始终存储1/3
a = 1/3; b = sym(1/3);a, b
例2-2 圆周率
求出
的300位有效数字。
vpa(pi, 300)
例2-3 试用符号型数据结构表示数值12345678901234567890
A = sym(12345678901234567890)
%% 显示结果:sym()内的数值被截断了
A =
12345678901234567168
应该这样:
B = sym('12345678901234567890')
%%结果正确
B =
12345678901234567890
02.02 矩阵与向量的输入方法
直接赋值语句: variable = expression
保留变量:ans,存放最近依次无赋值变量语言的运算结果
赋值语句的末尾加一个分号可以阻止显示运算结果。
例 2-5 矩阵输入方法
A = [1, 2, 3; 4, 5, 6; 7, 8, 9]
%% 动态定维
A = [1, 2, 3; 4 5, 6; 7, 8 9];
A = [[A; [1 2 3]], [1; 2; 3; 4]]
inv(A) % 求矩阵A的逆矩阵
例 2-6 复数矩阵输入
B = [1+9i 2+8i, 3+7i; 4+6i 5+5i, 6+4i; 7 + 3i , 8 + 2i, 0 + 1i]
% 注意空格,当1 +9i ,+前有空格时,会提示[串联的矩阵的维度不一致。]错误
B = [1 +9i 2+8i, 3+7i; 4+6i 5+5i, 6+4i; 7 + 3i , 8 + 2i, 0 + 1i]
% [1 + 9i]+前后空格正确
B = [1 + 9i 2+8i, 3+7i; 4+6i, 5+5i, 6+4i; 7 + 3i, 8 + 2i, 0 + 1i]
函数调用
函数调用语句
[returned.arguments] = function_name(input_arguments)
函数调用: [U S V] = svd(X)
函数可以通过不同方式被调用:
内核函数: *.m 函数
匿名函数、inline 函数(不建议使用)
重载函数、私有函数等
冒号表达式
默认步长为1
例2-7 冒号表达式生成
用不同的步距生成
间的向量
v1 = 0: 0.2: pi % 不包含pi
v1a = linspace(0, pi, 50) % 包含pi
v2 = 0: -0.1: pi % 空的 1×0 double 行矢量
v3 = 0: pi % 0 1 2 3
v4 = pi: -1: 0 %3.1416 2.1416 1.1416 0.1416
子矩阵提取
B = A(v1, v2)
v1 表示子矩阵要保留的行号构成的向量
v2表示子矩阵的列好构成的向量
:,表示要提取所有行或列,取决于其位置
end 的使用
% matlab 子矩阵获取
B = A(1:2:end, :), C = A([1 1 1 1], :)
02.03 矩阵的代数运算
矩阵表示: 矩阵A,n行m列,被称作n x m 矩阵
转置
% matlab 共轭转置(Hermite转置)
C = A'
% matlab 一般转置
C = A .'
加减法
数学表示:
% 加法matlab表示
C = A + B C = A - B
任一变量可为标量;如果矩阵A、B维度不同,matlab报错。
乘法
数学表示:
% 乘法 matlab 表示
C = A * B
除法
矩阵左除:
求解线性方程组:AX = B
% matlab 求解左除
X = A \B
最小二乘解;若A为非奇异方阵,则
使得误差最小:
矩阵右除:
求解线性方程组: XA = B
%矩阵右除
X = B / A
等效的运算: B / A = (A' \ B')'
矩阵翻转
左右翻转:
, B = fliplr(A)
上下翻转:
, B = flipud(A)
旋转90°:
, D = rot90(A)
如果旋转180°:D = rot90(A, k)
矩阵乘方
求矩阵A的x次幂:
数学描述:
% matlab 矩阵的乘方
F = A^x
A必须为方阵
x为整数
x为非整数
开方的多解:旋转
得出一个解,乘以k-1次标量
例2-9 矩阵的三次方根
% 矩阵的3次方根
A = [1, 2, 3; 4, 5, 6; 7, 8, 0];
C = A^(1/3), e = norm(A-C^3)
另两个根
j1 = exp(sqrt(-1)*2*pi/3);
A1 = C * j1, A2 = C * j1^2,
norm(A-A1^3), norm(A-A2^3)
点运算
矩阵对应元素的直接运算
例如:
C = A .* B,
B = A .^ A,
B = A .^ A =
A = [1, 2, 3; 4, 5, 6; 7, 8, 0]; B = A .^ A
绘图:
x = [-2, -1, 0, 1, 2], y = x .^2; plot(x, y)
02.04 矩阵的其它运算
逻辑运算 – 相应元素间的运算
与运算: A & B
或运算: A | B
非运算: B = ~A
异或运算: xor(A, B)
比较关系符:>, >=, <, <=, ==, ~=, find(), all(), any()
例2-10 比较运算实例
%比较运算实例
A = [1, 2, 3; 4, 5, 6; 7, 8, 0];
i = find(A >= 5)'
[i,j] = find(A >= 5)
a1 = all(A >= 5), a2 = any(A >= 5)
解析结果的化简与变换
函数simplify()用于数学公式的化简
= simplify(s)
其它常用化简函数
numden(), collect(), expand(), factor()
例2-11 多项式处理
化简多项式
用不同函数求解
% 多项式化简
syms s;
P = (s+3)^2 * (s^2 + 3*s + 2) * (s^3 + 12*s^2 + 48*s + 64)
P1 = simplify(P)
% 多项式因式分解
P3 = factor(P), P4 = prod(P3)
变量替换
转换成LaTex表示 – 需要LaTex环境解释
例2-12 双线性变换
试用 s = (z-1)/(z+1)对 P(s) 进行双线性变换
多项式表达式:
% matlab 变量替换代码
syms s z;
P = (s+3)^2*(s^2+3*s+2)*...
(s^3+12*s^2+48*s+64);
P1 = simplify(subs(P, s, (z-1)/(z+1)))
latex(P1)
% '\frac{8\, z\, {\left(2\, z + 1\right)}^2\, \left(3\, z + 1\right)\, {\left(5\, z + 3\right)}^3}{{\left(z + 1\right)}^7}'
基本数论运算
Function | Calling format |
---|---|
floor() | n = floor(x) |
ceil() | n = ceil(x) |
round() | n = round(x) |
fix() | n = fix(x) |
rat() | [n, d] = rat(x) |
rem() | B = rem(A, C) |
gcd() | k = gcd(n, m) |
lcm() | k = lcm(n, m) |
factor() | factor(n) |
isprime() | v1 = isprime(v) |
例 2-13 不同的取整函数
运用各种函数,对下面的数据进行取整运算
-0.2765, 0.5772, 1.4597, 2.1091, 1.191, -1.6187
% 各种取整
A = [-0.2765, 0.5772, 1.4597, 2.1091, 1.191, -1.6187];
v1 = floor(A), v2 = ceil(A),
v3 = round(A), v4 = fix(A)
例2-14 Hilbert矩阵
假设3x3的Hilbert矩阵可以右hilb()定义,试对其进行有理数变换。
% Hilbert矩阵有理数变换
A = hilb(3); [n, d] = rat(A)
例2-15 最大公约数和最小公倍数
给定两个整数1856120和1483720。
求其最大公约数和最小公倍数。
求其所得出的最小公倍数的质因数分解。
m = sym(1856120); n = sym(1483720);
gcd(m, n), lcm(m, n),
factor(lcm(n, m))
例2-16 找出1000以内全部质数
% 找出1000以内所有质数
A = 1:1000; B = A(isprime(A))
例2-17 全排列计算
全排列函数: perms()
给出1到5的全排列.
% matlab 全排列
P = perms(1:5), size(P)
P = perms('abcde'), size(P)
02.05 matlab 流程结构
循环语句
for循环 : for i = v, loop programs body, end
while循环:
while (condition)
loop structure body,
end
例2-18 求解
% for loop
s1 = 0; for i = 1: 100, s1 = s1 + i; end, s1
% while loop
s2 = 0; i = 1;
while(i <= 100), s2 = s2 + i; i = i + 1; end, s2
% 简单sum函数
sum(1:100)
例2-19 while循环
用循环求解最小的m,使下式成立
s = 0; m = 0;
while (s <= 10000), m = m + 1;
s = s + m; end, s, m
例2-20 向量化编程
求和
% 循环求解
tic, s = 0;
for i = 1: 100000, s = s + 1/2^i + 1/3^i;
end; toc
% 时间已过 0.038848 秒。
% %%%%%%%%%%%
% 向量化编程
tic, i = 1:100000;
s = sum(1 ./2 .^i + 1./ 3 .^i); toc
% 时间已过 0.019701 秒。
条件语句
if (condition 1)
statement group 1
elseif (condition 2)
statement group 2
else
statement group n + 1
end
例2-21 for循环+if
用for循环求解最大的m,试下式成立:
s = 0;
for i = 1: 10000, s = s + i;
if s > 10000, break;
end,
end, i, s
开关语句
switch switch expression
case expression 1, statements 1
case {expression2, expression 3, ..., expression m}, statements 2
...
otherwise, statements n
end
try catch
try,
statement group1,
catch,
statement group2,
end
02.06 函数编写
例2-22 M-脚本文件实现
问题:
m-脚本以m-文件形式存取
% 以下内容保持到文件s10000.m
s = 0; m = 0;
while (s <= 10000),
m = m + 1,
s = s + m;
end,
s, m
% 直接执行s10000,即调用s10000.m脚本运行
% 函数基本结构
function [return vars] = funcname(input vars)
comments led by %
input and output variables check
main body of the function
例2-23 m-函数实现
% 函数实例
function [ m, s ] = findsum( k )
%findsum summary
% findsum Detailed description
s = 0; m = 0;
while (s <= k)
m = m + 1;
s = s + m;
end;
end
% call function
[m1, s1] = findsum(145323)
例2-24 Hilbert长方形矩阵
编写一个函数生成 n x m Hilbert 矩阵
function [ output_args ] = myhilb( n, m )
if nargin == 1
m = n;
end;
for i = 1:n,
for j = 1: m
H(i, j) = 1/(i+j-1);
end,
end;
end
% call function
% H = myhilb(sym(4), 3)
% H = myhilb(3)
例2-25 阶乘的递归调用 n! = n(n-1)!
% n!
function k = my_fact(n)
if nargin ~= 1
error('Error: Only one input variable accepted');
end;
if abs(n - floor(n)) > eps || n < 0
error('n should be a non-negative integer');
end;
if n > 1, k = n * my_fact(n-1);
elseif any([0, 1] == n)
k =1;
end;
end
阶乘的不同计算方法
factorial(n)
prod(1:n)
gamma(1+n)
gamma(1 + sym(n))
可变输入输出
输入输出变量: varargin, varargout
变量的提取 - 单元数组(cell)
varargin{1}, varargin{2}, …, varargin{n}
例2-27 任意多输入变量
conv()可计算两个多项式的积,试用varargin实现任意多个多项式的积
function a = convs(varargin), a = 1;
for i = 1: nargin;
a = conv(a, varargin{i});
end;
end
% call convs
% P = [1 2 4 0 5]; Q = [1 2];
% F = [1 2 3]; D = convs(P, Q, F), E = conv(conv(P, Q), F)
% G = convs(P, Q, F, [1, 1], [1, 3], [1, 1])
匿名函数
f = @(list of variables) function_contents
伪代码与代码保密
伪代码:提高程序的执行速度,保密,把ascii的.m文件转换成二进制文件
伪代码语句:
pcode mytest
pcode *.m
pcode mytest -inplace
一定要先保存好.m源文件,pcode没有对应恢复命令。没办法通过.p文件恢复成.m文件。
02.07 二维曲线的绘制
例2-28 函数曲线绘制与检验
绘制函数:
x = [-pi: 0.05: pi];
y = sin(tan(x)) - tan(sin(x)); plot(x, y)
不同步距效果
x = [-pi: 0.0001: pi];
y = sin(tan(x)) - tan(sin(x)); plot(x, y)
例2-29 分段函数
绘制饱和函数方程
互斥条件:
x = [-2: 0.02: 2];
y = 1.1 * sign(x) .* (abs(x) > 1.1) + ...
x .* (abs(x) <= 1.1); plot(x, y)
更简单的命令–折线
plot([-2, -1.1, 1.1, 2], [-1.1, -1.1, 1.1, 1.1])
plot其它调用格式
t仍为向量,而y为矩阵,
plot(t, y)
t和y为矩阵,且t和y的行和列数均相同
假设有多对这样的向量或矩阵
plot
更一般的调用格式
改变曲线性质
h = plot(
, option 1,
, option 2,
, option m)
图形修饰与属性设置
每个窗口、曲线、坐标轴都是一个对象
对象的属性可以通过set()来设置
可以通过函数get()来获取
set(handle, ‘p_name1’, …, p_value1, ‘p_name2’, p_value2, …)
v = get (object, ‘p_name’)
图形对象的属性可以通过快捷菜单(鼠标右键)直接修改
多纵轴曲线的绘制
例2-30 试绘制曲线
x = 0: 0.01: 2*pi;
y1 = sin(x); y2 = 0.01 * cos(x);
plot(x, y1, x, y2, '--')
x = 0: 0.01: 2*pi;
y1 = sin(x); y2 = 0.01 * cos(x);
plotyy(x, y1, x, y2)
三、四纵轴图形可以下载相应函数绘制
plotyyy(), plot4y(), 从MathWorks File Exchange 下载,还可以使用 plotxx()函数
02.08 特殊二维图形
其它二维图形绘制语言
* | * |
---|---|
bar(x, y) | comet(x, y) |
compass(x,y) | errorbar(x, y, ) |
feather(x, y) | fill(x, y, c) |
hist(y, n) | loglog(x, y) |
polar(x, y) | quiver(x, y) |
stairs(x, y) | stem(x, y) |
semilogx(x, y) | semilogy(x,y) |
例2-31 极坐标曲线
绘制极坐标函数:
theta = 0: 0.01: 6*pi;
rho = 5*sin(4*theta/3);
polar(theta, rho)
rho = 5*sin(theta/3);
polar(theta, rho)
例2-32 不同的二维图形
t = 0: .2: 2*pi; y = sin(t);
subplot(2, 2, 1), stairs(t, y)
subplot(2, 2, 2), stem(t, y)
subplot(2, 2, 3), bar(t, y)
subplot(2, 2, 4), semilogx(t, y)
绘制区域设置
分割不同区域: subplot(n, m, k)
隐函数绘制及应用
隐函数:
隐函数绘图语句: ezplot(implicit function expression)
用字符串表示隐函数
默认区域是
其它语法: ezplot(im_function, [
])
ezplot(im_function, [
])
例2-33 隐函数曲线
试绘制隐函数:
% 默认范围 -2*pi, 2*pi
h = ezplot('x^2*sin(x+y^2) + y^2*exp(x+y)+5*cos(x^2+y)');
set(h, 'Color', 'b')
% 扩大范围
h = ezplot('x^2*sin(x+y^2) + y^2*exp(x+y) + 5*cos(x^2+y)', [-10, 10]);
set(h, 'Color', 'r');
图形修饰
直接采用工具栏
文字修饰
特殊符号表 LaTex
上、下标分别用 ^ 和 _ 表示
a_2^2 + b_2^2 = c_2^2 ?
数据文件的读取与存储
可以采用save和load命令存储和读取数据
save mydat A B C
save /ascii mydat.dat A B C
X = load(filename)
matlab和excel交互数据
X = xlsread(filename, range) ‘B6:C67’
写文件: xlswrite()
例2-34 Excel文件读取
已知Excel文件 census.xls 给出某省人口数
第5-67行存储数据
B列存储年份,C列存储人口数
先读入matlab再绘图
X = xlsread('census.xls', 'B5:C67');
t = X(:1); p = X(:, 2);
plot(t, p)
更加单方法 – Copy & Paste
02.09 三维图形绘制
三维曲线绘制
三维曲线绘制(空间质点的运动轨迹)
plot3(x, y, z)
plot3(
, option 1,
, option 2,
option m)
其它三维曲线绘制函数
stem3, fill3, bar3等
绘制三维轨迹图: comet3
例2-35 空间质点的位置
试绘制参数方程:
质点的空间位置(随时间变化)
其中,
t = 0: 0.01: 2*pi;
x = t .^ 3 .* sin(3*t) .* exp(-t);
y = t .^ 3 .* cos(3*t) .* exp(-t);
z = t .^ 2;
plot3(x, y, z), grid
其它曲线绘制
% 其它曲线绘制
stem3(x, y, z); hold on;
plot3(x, y, z); grid; hold off
% 动画效果
figure; comet3(x, y, z)
图形窗口的工具栏
3D绘图和视角变换
读取坐标值、局部放大
三维曲面绘制
一般曲面绘制: z = f(x, y)
[x, y] = meshgrid(
)
z = …, for instance z = x .* y
surf(x, y, z) or mesh(x, y, z)
其它函数surfl(), surfc()
等高线绘制
contour(), contours()
例2-36 三维曲面
给出二元函数如下, 绘制3D图形
[x, y] = meshgrid(-3: 0.1: 3, -2: 0.1: 2);
z = (x .^ 2 - 2*x) .* exp(-x .^ 2 -y .^ 2 - x .* y);
mesh(x, y, z)
% 表面图:
surf(x, y, z)
例2-37 试绘制出二元函数
绘制3D图形
[x, y] = meshgrid(-2: .1: 2);
z = 1 ./ (sqrt((1-x).^2 + y .^2)) + 1 ./ (sqrt((1+x).^2 + y .^2));
surf(x, y, z), shading flat
例2-38 分段函数处理
绘制如下分段函数的三维曲面
分段函数求值: 互斥的不等式、点运算
[x1, x2] = meshgrid(-1.5: .1: 1.5, -2: .1: 2);
p = 0.5457*exp(-0.75*x2.^2 - 3.75*x1.^2-1.5*x1).*(x1+x2>1) + ...
0.7575*exp(-x2.^2-6*x1.^2).*((x1+x2>-1)&(x1+x2<=1)) +...
0.5457*exp(-0.75*x2.^2-3.75*x1.^2+1.5*x1).*(x1+x2<=-1);
surf(x1, x2, p), xlim([-1.5 1.5]);
参数方程的表面图
三维函数参数方程
参数区间:
ezsurf(
)
默认区间:
例2-40 Mobius带的绘制
数学模型:
% Mobius带
syms u v; x= cos(u) + v*cos(u)*cos(u/2);
y = sin(u) + v*sin(u)*cos(u/2);
z = v*sin(u/2);
ezsurf(x, y, z, [0, 2*pi, -0.5, 0.5])
球面绘制
matlab函数生成数据[x, y, z] = shpere(n)
生成矩阵(n+1) x (n+1)
例2-41 球面绘制
绘制两个球, 一个圆心在原点, 半径为1,另一个圆心(0.9, -0.8, 0.6),半径为0.3
[x, y, z] = sphere(50); surf(x, y, z), hold on
x1 = 0.3*x + 0.9; y1 = 0.3*y - 0.8; z1 = 0.3*z + 0.6;
surf(x1, y1, z1), hold off
柱面绘制
某曲线r沿z轴旋转一周得出柱面
[x, y, z] = cylinder(r, n)
例2-42 柱面曲线方程
z0 = -1: 0.1: 3; r = exp(-z0.^2/2).*sin(z0);
[x, y, z] = cylinder(r); z = -1 + 4*z; surf(x, y, z)
02.10 特殊三维图形
等高线绘制
各种等高线绘制方法
直接绘制: contour(x, y, z, n)
带有标注的等高线:
[C, h] = countour(x, y, z, n)
clabel(C, h)
三维等高线countour3(x, y, z, n) coutourf(x, y, z, n)
例2-43 等高线
例2-38的分段函数,绘制各种等高线
生成数据,再绘制图形
[x1, x2] = meshgrid(-1.5: .1: 1.5, -2: .1: 2);
p = 0.5457*exp(-0.75*x2.^2 - 3.75*x1.^2-1.5*x1).*(x1+x2>1) + ...
0.7575*exp(-x2.^2-6*x1.^2).*((x1+x2>-1)&(x1+x2<=1)) +...
0.5457*exp(-0.75*x2.^2-3.75*x1.^2+1.5*x1).*(x1+x2<=-1);
[C, h] = contour(x1, x2, p); clabel(C, h)
contourf(x1, x2, p);
figure; contour3(x1, x2, p, 30)
三维隐函数图绘制
三元隐函数:
下载 ezimplot3()函数
ezimplot3(fun, [
])
默认范围:
fun可以为[隐函数字符串、匿名函数、M-函数]
例2-44 三元隐函数绘制
已知三元隐函数
f = 'x*sin(y+z^2) + y^2*cos(x+z)+z*x*cos(z+y^2)';
f = @(x, y, z)x*sin(y+z^2) + y^2*cos(x+z)+z*x*cos(z+y^2);
ezimplot3(f, [-1, 1]);
f1 = 'x^2+y^2+z^2-1'; ezimplot3(f1, [-1 1]);
三维图形视角设置
两种方法改变图形视角
直接采用工具栏
命令语句:view()
读角度:
定义为方位角,
定义为仰角
例 2-45 三视图的设置
函数:
默认视角的提取: [a, b] = view(3);
三视图绘制:
[x, y] = meshgrid(-3: 0.1: 3, -2: 0.1: 2);
z = (x .^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
subplot(224), surf(x,y,z),
subplot(221), surf(x,y,z), view(0, 90);
subplot(222), surf(x,y,z), view(90, 0);
subplot(223), surf(x,y,z), view(0, 0);
三维曲面旋转
三维曲面旋转函数: rotate(h, v,
)
其中: 三维曲线句柄,旋转基线,旋转角度
如: 绕x轴正向旋转, v = [1, 0, 0]
例2-46 三维图形旋转
分段函数曲面旋转绘制
[x, y] = meshgrid(-1: 0.04: 1, -2: 0.04: 2);
z = 0.5457*exp(-0.75*y.^2 - 3.75*x.^2-1.5*x).*(x+y>1) + ...
0.7575*exp(-y.^2-6*x.^2).*((x+y>-1)&(x+y<=1)) +...
0.5457*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);
h = surf(x, y, z);
% 绕x轴正方向旋转,基线 v = [1, 0, 0]
rot_ax = [1, 0, 0]; rotate(h, rot_ax, 15)
% 绕基线 v = [1, 1, 1] 旋转
figure;
h = surf(x, y, z); rot_ax = [1, 1, 1];
rotate(h, rot_ax, 15)
旋转的动画演示
绕x轴正向旋转360度的动画
每步(0.01秒)旋转1度
循环结构旋转
figure;h = surf(x,y,z);
r_ax = [1, 0, 0]; axis tight
for i = 0:360
rotate(h, r_ax, 1); pause(0.01);
end;
02.11 四维图形绘制
四维图形
时间驱动的三维图形的动画
三维图形的切面 – 体视化(volume visualization)
体视化举例: CT、固体内部温度、密度,需要用切面观察
流体的流速、液体的浓度分布
matlab函数:slice(x, y, z, V,
)
例2-47 三维实体的切片
三元函数
普通切面
生成网格数据
计算网格各点的函数值(点运算)
绘图
[x, y, z] = meshgrid(0: 0.1: 2);
V = sqrt(x .^ x + y .^((x+y)/2) + z .^((x+y+z)/3) );
slice(x, y, z, V, [1 2], [1 2], [0 1]);
特殊切面
先构造 z = 1 平面
将该平面沿x轴旋转45度
用该切面给原体视化图形切片
[x0, y0] = meshgrid(0: 0.1: 2);
z0 = ones(size(x0));
h = surf(x0, y0, z0);
rotate(h, [1,0,0], 45);
x1 = get(h, 'XData');
y1 = get(h, 'YData');
z1 = get(h, 'ZData');
slice(x, y, z, V, x1, y1, z1),
hold on, slice(x, y, z, V, 2, 2, 0), hold off
shading interp
编写体视化程序界面
程序界面 vol_visual4d()
数学函数: V = f(x, y, z)
vol_visual.fig, vol_visual4d.m
调用方法
生成体视化数据:x, y, z, V
调用该函数: vol_visual4d(x, y, z, V)
通过界面任意设置切面
例2-48 利用界面体视化
三元函数
体视化切面绘图
通过滚动杆调整切面位置
通过捡取框选择是否显示某轴切面
[x,y,z] = meshgrid(0: 0.1: 2);
V = sqrt(x .^ x + y .^ ((x+y)/2) + z .^ ((x+y+z)/3) );
vol_visual4d(x, y, z, V);