Tag Archives: matlab
how can i convert 2d plot in to 3d surface ??
I have this small code that plots ogive shape nose in 2d.
I want to develop 3d surface from this code !!! please help and guide….
I am using surf command and tried to turn vector in to matrix….
L = 500;
R = 79;
n= 0.4;
x = 0:0.01:50;
y= R*(x/L).^n;
z= 1.*(x/L).^n;
plot(x,y,’r’,’linewidth’,2)
hold on
plot(x,-y,’r’,’linewidth’,2)
xlabel(‘Length’)
ylabel(‘Radius’)
title(‘Nose Profile’)
axis(‘equal’);I have this small code that plots ogive shape nose in 2d.
I want to develop 3d surface from this code !!! please help and guide….
I am using surf command and tried to turn vector in to matrix….
L = 500;
R = 79;
n= 0.4;
x = 0:0.01:50;
y= R*(x/L).^n;
z= 1.*(x/L).^n;
plot(x,y,’r’,’linewidth’,2)
hold on
plot(x,-y,’r’,’linewidth’,2)
xlabel(‘Length’)
ylabel(‘Radius’)
title(‘Nose Profile’)
axis(‘equal’); I have this small code that plots ogive shape nose in 2d.
I want to develop 3d surface from this code !!! please help and guide….
I am using surf command and tried to turn vector in to matrix….
L = 500;
R = 79;
n= 0.4;
x = 0:0.01:50;
y= R*(x/L).^n;
z= 1.*(x/L).^n;
plot(x,y,’r’,’linewidth’,2)
hold on
plot(x,-y,’r’,’linewidth’,2)
xlabel(‘Length’)
ylabel(‘Radius’)
title(‘Nose Profile’)
axis(‘equal’); surf, 2d, plotting, 3d plots MATLAB Answers — New Questions
Sorting a Cell Array
Hi,
Sorry if this has been asked / answered many times or very simple etc. etc!
If I have a cell array say Array = {B C A;1 3 2; D E F}
Is it possible to sort the first row in order, so its A B C, with the second and third rows changing corresponding to to this initial sort?
Final answer being {A B C;2 1 3;F D E}
Thanks in advance
MatthewHi,
Sorry if this has been asked / answered many times or very simple etc. etc!
If I have a cell array say Array = {B C A;1 3 2; D E F}
Is it possible to sort the first row in order, so its A B C, with the second and third rows changing corresponding to to this initial sort?
Final answer being {A B C;2 1 3;F D E}
Thanks in advance
Matthew Hi,
Sorry if this has been asked / answered many times or very simple etc. etc!
If I have a cell array say Array = {B C A;1 3 2; D E F}
Is it possible to sort the first row in order, so its A B C, with the second and third rows changing corresponding to to this initial sort?
Final answer being {A B C;2 1 3;F D E}
Thanks in advance
Matthew array, sort MATLAB Answers — New Questions
Index exceeds the number of array elements. Index must not exceed 3. HELP
model code:
clear all;
close all;
S = 99;
I = 1;
R = 0;
N = 100; %Total population
beta= 0.1; % birth rate
alpha= 0.1; % infection person to person rate
lambda= 0.3; % infection by water rate
vac= 0.05; % recovery by vaccination rate
d= 0.03; % death rate
gamma= 0.8; % recovery rate
c= 0.9; % rate of contamination
m= 0.4; % rate of decay of V. cholerae
B= 0.0; % initial concentration of V. cholerae
t_f = 500; %Ending time of simulation
Q = [beta alpha lambda vac d gamma c m B];
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q);
figure(1)
plot(t,y(:,1),’k-‘,t,y(:,2),’r–‘,t,y(:,3),’b:’);
xlabel(‘bf Time (days)’);
ylabel(‘bf Number of People by Category’);
legend(‘S’,’I’,’R’);
z=y(end,:)’
SN=y(:,1)/N;
IN=y(:,2)/N;
figure(2)
plot(IN,SN);
xlabel(‘bf I/N’);
ylabel(‘bf S/N’);
r0=beta/beta
[maxIN,y_maxtime]=max(y(:,2)/N);
maxIN
maxtime=y(y_maxtime)
eqIN=y(100,2)/N;
eqIN
function code:
function dy=cholera_de(t,Y,flag,Q)
beta= Q(1);
alpha= Q(2);
lambda= Q(3);
vac= Q(4);
d= Q(5);
gamma= Q(6);
c= Q(7);
m= Q(8);
S= Y(1);
I= Y(2);
R= Y(3);
B=Y(4);
N= S+I+R;
dy(1,1)= beta – alpha*I – lambda*B – vac*S – d;
dy(2,1)= alpha*I + lambda*B – d – gamma*I;
dy(3,1)= gamma*I + vac*S – d;
dy(4,1)= c*I – m*B;
Gives me this error:
Index exceeds the number of array elements. Index must not exceed 3.
Error in cholera_de (line 15)
B=Y(4);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in cholera_model (line 23)
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q);model code:
clear all;
close all;
S = 99;
I = 1;
R = 0;
N = 100; %Total population
beta= 0.1; % birth rate
alpha= 0.1; % infection person to person rate
lambda= 0.3; % infection by water rate
vac= 0.05; % recovery by vaccination rate
d= 0.03; % death rate
gamma= 0.8; % recovery rate
c= 0.9; % rate of contamination
m= 0.4; % rate of decay of V. cholerae
B= 0.0; % initial concentration of V. cholerae
t_f = 500; %Ending time of simulation
Q = [beta alpha lambda vac d gamma c m B];
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q);
figure(1)
plot(t,y(:,1),’k-‘,t,y(:,2),’r–‘,t,y(:,3),’b:’);
xlabel(‘bf Time (days)’);
ylabel(‘bf Number of People by Category’);
legend(‘S’,’I’,’R’);
z=y(end,:)’
SN=y(:,1)/N;
IN=y(:,2)/N;
figure(2)
plot(IN,SN);
xlabel(‘bf I/N’);
ylabel(‘bf S/N’);
r0=beta/beta
[maxIN,y_maxtime]=max(y(:,2)/N);
maxIN
maxtime=y(y_maxtime)
eqIN=y(100,2)/N;
eqIN
function code:
function dy=cholera_de(t,Y,flag,Q)
beta= Q(1);
alpha= Q(2);
lambda= Q(3);
vac= Q(4);
d= Q(5);
gamma= Q(6);
c= Q(7);
m= Q(8);
S= Y(1);
I= Y(2);
R= Y(3);
B=Y(4);
N= S+I+R;
dy(1,1)= beta – alpha*I – lambda*B – vac*S – d;
dy(2,1)= alpha*I + lambda*B – d – gamma*I;
dy(3,1)= gamma*I + vac*S – d;
dy(4,1)= c*I – m*B;
Gives me this error:
Index exceeds the number of array elements. Index must not exceed 3.
Error in cholera_de (line 15)
B=Y(4);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in cholera_model (line 23)
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q); model code:
clear all;
close all;
S = 99;
I = 1;
R = 0;
N = 100; %Total population
beta= 0.1; % birth rate
alpha= 0.1; % infection person to person rate
lambda= 0.3; % infection by water rate
vac= 0.05; % recovery by vaccination rate
d= 0.03; % death rate
gamma= 0.8; % recovery rate
c= 0.9; % rate of contamination
m= 0.4; % rate of decay of V. cholerae
B= 0.0; % initial concentration of V. cholerae
t_f = 500; %Ending time of simulation
Q = [beta alpha lambda vac d gamma c m B];
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q);
figure(1)
plot(t,y(:,1),’k-‘,t,y(:,2),’r–‘,t,y(:,3),’b:’);
xlabel(‘bf Time (days)’);
ylabel(‘bf Number of People by Category’);
legend(‘S’,’I’,’R’);
z=y(end,:)’
SN=y(:,1)/N;
IN=y(:,2)/N;
figure(2)
plot(IN,SN);
xlabel(‘bf I/N’);
ylabel(‘bf S/N’);
r0=beta/beta
[maxIN,y_maxtime]=max(y(:,2)/N);
maxIN
maxtime=y(y_maxtime)
eqIN=y(100,2)/N;
eqIN
function code:
function dy=cholera_de(t,Y,flag,Q)
beta= Q(1);
alpha= Q(2);
lambda= Q(3);
vac= Q(4);
d= Q(5);
gamma= Q(6);
c= Q(7);
m= Q(8);
S= Y(1);
I= Y(2);
R= Y(3);
B=Y(4);
N= S+I+R;
dy(1,1)= beta – alpha*I – lambda*B – vac*S – d;
dy(2,1)= alpha*I + lambda*B – d – gamma*I;
dy(3,1)= gamma*I + vac*S – d;
dy(4,1)= c*I – m*B;
Gives me this error:
Index exceeds the number of array elements. Index must not exceed 3.
Error in cholera_de (line 15)
B=Y(4);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in cholera_model (line 23)
[t,y]=ode45(‘cholera_de’,[0:t_f/100:t_f],[S I R]’,[],Q); ode45, index MATLAB Answers — New Questions
Error Using Reshape in MATLAB
How do I correct the following error in MATLAB? MATLAB code attached.
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray), [], 1);
% BPSK modulation
Eb_No_low = 0; % Low SNR
Eb_No_high = 4; % High SNR
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(double(lenna_bits), SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_low); title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(double(lenna_bits), SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_high); title(‘Received Image (4 dB SNR)’);
% Linear error detection code
% Example: Hamming (7,4) code
parityMatrix = [1 1 1 0 1 0 0; 1 1 0 1 0 1 0; 1 0 1 1 0 0 1];
generatorMatrix = [eye(4) parityMatrix’];
% Encode the data
encoded_data = mod(lenna_bits * generatorMatrix, 2);
% Add noise for linear error detection code
received_data = awgn(double(encoded_data), Eb_No_low, ‘measured’);
% Syndrome lookup table for error detection
syndrome_table = syndtable(parityMatrix);
% Decoding with error detection
decoded_data = zeros(size(encoded_data));
errors = zeros(size(encoded_data, 1), 1);
for i = 1:size(encoded_data, 1)
syndrome = mod(received_data(i, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(i) = 1;
else
decoded_data(i, 🙂 = received_data(i, :);
end
end
% Count number of retransmission requests at different SNRs
SNRs = [0, 2, 4, 6, 8, 10];
retransmissions = zeros(size(SNRs));
for i = 1:length(SNRs)
SNR = 10^(SNRs(i)/10);
received_data = awgn(encoded_data, SNR, ‘measured’);
errors = zeros(size(encoded_data, 1), 1);
for j = 1:size(encoded_data, 1)
syndrome = mod(received_data(j, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(j) = 1;
retransmissions(i) = retransmissions(i) + 1;
end
end
end
% Plot number of retransmissions against SNR values
figure;
plot(SNRs, retransmissions, ‘-o’);
xlabel(‘SNR (dB)’);
ylabel(‘Number of Retransmissions’);
title(‘Number of Retransmissions vs SNR’);
% Error correction code
% Example: Reed-Solomon code
n = 255;
k = 223;
t = 16;
rs_encoder = comm.RSEncoder(n, k);
rs_decoder = comm.RSDecoder(n, k);
% Encode data
encoded_rs = step(rs_encoder, double(lenna_bits));
% Add noise for Reed-Solomon code
received_rs_low = awgn(double(encoded_rs), Eb_No_low, ‘measured’);
received_rs_high = awgn(double(encoded_rs), Eb_No_high, ‘measured’);
% Decode received data
decoded_rs_low = step(rs_decoder, received_rs_low);
decoded_rs_high = step(rs_decoder, received_rs_high);
% Reshape decoded bits to original image size
decoded_image_rs_low = reshape(decoded_rs_low, size(lenna_gray));
decoded_image_rs_high = reshape(decoded_rs_high, size(lenna_gray));
% Plot original and received images with error correction
figure;
subplot(1,2,1); imshow(decoded_image_low); title(‘Received Image (0 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_low); title(‘Received Image (0 dB SNR, Error Correction)’);
figure;
subplot(1,2,1); imshow(decoded_image_high); title(‘Received Image (4 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_high); title(‘Received Image (4 dB SNR, Error Correction)’);
Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that
dimension.
Error in Computerassignmentfinal (line 26)
decoded_image_low = reshape(decoded_low, size(lenna_gray));How do I correct the following error in MATLAB? MATLAB code attached.
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray), [], 1);
% BPSK modulation
Eb_No_low = 0; % Low SNR
Eb_No_high = 4; % High SNR
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(double(lenna_bits), SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_low); title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(double(lenna_bits), SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_high); title(‘Received Image (4 dB SNR)’);
% Linear error detection code
% Example: Hamming (7,4) code
parityMatrix = [1 1 1 0 1 0 0; 1 1 0 1 0 1 0; 1 0 1 1 0 0 1];
generatorMatrix = [eye(4) parityMatrix’];
% Encode the data
encoded_data = mod(lenna_bits * generatorMatrix, 2);
% Add noise for linear error detection code
received_data = awgn(double(encoded_data), Eb_No_low, ‘measured’);
% Syndrome lookup table for error detection
syndrome_table = syndtable(parityMatrix);
% Decoding with error detection
decoded_data = zeros(size(encoded_data));
errors = zeros(size(encoded_data, 1), 1);
for i = 1:size(encoded_data, 1)
syndrome = mod(received_data(i, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(i) = 1;
else
decoded_data(i, 🙂 = received_data(i, :);
end
end
% Count number of retransmission requests at different SNRs
SNRs = [0, 2, 4, 6, 8, 10];
retransmissions = zeros(size(SNRs));
for i = 1:length(SNRs)
SNR = 10^(SNRs(i)/10);
received_data = awgn(encoded_data, SNR, ‘measured’);
errors = zeros(size(encoded_data, 1), 1);
for j = 1:size(encoded_data, 1)
syndrome = mod(received_data(j, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(j) = 1;
retransmissions(i) = retransmissions(i) + 1;
end
end
end
% Plot number of retransmissions against SNR values
figure;
plot(SNRs, retransmissions, ‘-o’);
xlabel(‘SNR (dB)’);
ylabel(‘Number of Retransmissions’);
title(‘Number of Retransmissions vs SNR’);
% Error correction code
% Example: Reed-Solomon code
n = 255;
k = 223;
t = 16;
rs_encoder = comm.RSEncoder(n, k);
rs_decoder = comm.RSDecoder(n, k);
% Encode data
encoded_rs = step(rs_encoder, double(lenna_bits));
% Add noise for Reed-Solomon code
received_rs_low = awgn(double(encoded_rs), Eb_No_low, ‘measured’);
received_rs_high = awgn(double(encoded_rs), Eb_No_high, ‘measured’);
% Decode received data
decoded_rs_low = step(rs_decoder, received_rs_low);
decoded_rs_high = step(rs_decoder, received_rs_high);
% Reshape decoded bits to original image size
decoded_image_rs_low = reshape(decoded_rs_low, size(lenna_gray));
decoded_image_rs_high = reshape(decoded_rs_high, size(lenna_gray));
% Plot original and received images with error correction
figure;
subplot(1,2,1); imshow(decoded_image_low); title(‘Received Image (0 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_low); title(‘Received Image (0 dB SNR, Error Correction)’);
figure;
subplot(1,2,1); imshow(decoded_image_high); title(‘Received Image (4 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_high); title(‘Received Image (4 dB SNR, Error Correction)’);
Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that
dimension.
Error in Computerassignmentfinal (line 26)
decoded_image_low = reshape(decoded_low, size(lenna_gray)); How do I correct the following error in MATLAB? MATLAB code attached.
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray), [], 1);
% BPSK modulation
Eb_No_low = 0; % Low SNR
Eb_No_high = 4; % High SNR
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(double(lenna_bits), SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_low); title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(double(lenna_bits), SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1); imshow(lenna_gray); title(‘Original Image’);
subplot(1,2,2); imshow(decoded_image_high); title(‘Received Image (4 dB SNR)’);
% Linear error detection code
% Example: Hamming (7,4) code
parityMatrix = [1 1 1 0 1 0 0; 1 1 0 1 0 1 0; 1 0 1 1 0 0 1];
generatorMatrix = [eye(4) parityMatrix’];
% Encode the data
encoded_data = mod(lenna_bits * generatorMatrix, 2);
% Add noise for linear error detection code
received_data = awgn(double(encoded_data), Eb_No_low, ‘measured’);
% Syndrome lookup table for error detection
syndrome_table = syndtable(parityMatrix);
% Decoding with error detection
decoded_data = zeros(size(encoded_data));
errors = zeros(size(encoded_data, 1), 1);
for i = 1:size(encoded_data, 1)
syndrome = mod(received_data(i, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(i) = 1;
else
decoded_data(i, 🙂 = received_data(i, :);
end
end
% Count number of retransmission requests at different SNRs
SNRs = [0, 2, 4, 6, 8, 10];
retransmissions = zeros(size(SNRs));
for i = 1:length(SNRs)
SNR = 10^(SNRs(i)/10);
received_data = awgn(encoded_data, SNR, ‘measured’);
errors = zeros(size(encoded_data, 1), 1);
for j = 1:size(encoded_data, 1)
syndrome = mod(received_data(j, 🙂 * parityMatrix’, 2);
if sum(syndrome) ~= 0 % Error detected
errors(j) = 1;
retransmissions(i) = retransmissions(i) + 1;
end
end
end
% Plot number of retransmissions against SNR values
figure;
plot(SNRs, retransmissions, ‘-o’);
xlabel(‘SNR (dB)’);
ylabel(‘Number of Retransmissions’);
title(‘Number of Retransmissions vs SNR’);
% Error correction code
% Example: Reed-Solomon code
n = 255;
k = 223;
t = 16;
rs_encoder = comm.RSEncoder(n, k);
rs_decoder = comm.RSDecoder(n, k);
% Encode data
encoded_rs = step(rs_encoder, double(lenna_bits));
% Add noise for Reed-Solomon code
received_rs_low = awgn(double(encoded_rs), Eb_No_low, ‘measured’);
received_rs_high = awgn(double(encoded_rs), Eb_No_high, ‘measured’);
% Decode received data
decoded_rs_low = step(rs_decoder, received_rs_low);
decoded_rs_high = step(rs_decoder, received_rs_high);
% Reshape decoded bits to original image size
decoded_image_rs_low = reshape(decoded_rs_low, size(lenna_gray));
decoded_image_rs_high = reshape(decoded_rs_high, size(lenna_gray));
% Plot original and received images with error correction
figure;
subplot(1,2,1); imshow(decoded_image_low); title(‘Received Image (0 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_low); title(‘Received Image (0 dB SNR, Error Correction)’);
figure;
subplot(1,2,1); imshow(decoded_image_high); title(‘Received Image (4 dB SNR, No Error Correction)’);
subplot(1,2,2); imshow(decoded_image_rs_high); title(‘Received Image (4 dB SNR, Error Correction)’);
Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that
dimension.
Error in Computerassignmentfinal (line 26)
decoded_image_low = reshape(decoded_low, size(lenna_gray)); reshape error MATLAB Answers — New Questions
Compute Confidence Interval about q of logistic binary fit
I have results from a binary study. The study was analyzied by others using Gonogo (https://www2.isye.gatech.edu/~jeffwu/sensitivitytesting/), the prefered code developed by/for the government for such studies. My code for anlaysis is shown below. I am able to fit a logistic curve using glmfit (magentia line in plot). I then compute the confidence limits about the probablility p using glmval (black lines). These black line match well with the results given by Gonogo (dashed blue lines). However, I am unable to figure out how to compute the confidence limits about the quantile q (the Gonongo results are shown as dashed red lines). I tried using coefCI, with the results in cyan. It is doing something but I don’t understand exactly what, and it does not really match the dashed red lines that I am looking for. Does anyone know how I can compute the confidence limits about the quantile q? Your assistance is appreciated.
%% Output from Gonogo – Main code in next section
R_mu = 10.1707909;
R_sigma = 0.9344091;
R_b0 = -10.88473;
R_b1 = 1.070195;
R_output = [
1.124021, 5.5, 9.875979, 0, 0, 0.3761890
1.321891, 5.611111, 9.900331, 0, 1.e-06, 0.3861210
1.519710, 5.722222, 9.924735, 0, 1.e-06, 0.3961490
1.717474, 5.833333, 9.949193, 0, 2.e-06, 0.4062690
1.915181, 5.944444, 9.973708, 0, 3.e-06, 0.4164760
2.112825, 6.055556, 9.998286, 0, 5.e-06, 0.4267660
2.310404, 6.166667, 10.022929, 0, 9.e-06, 0.4371330
2.507913, 6.277778, 10.047643, 0, 1.5e-05, 0.4475740
2.705346, 6.388889, 10.072432, 0, 2.6e-05, 0.4580840
2.902698, 6.5, 10.097302, 0, 4.3e-05, 0.4686570
3.099963, 6.611111, 10.122259, 0, 7.e-05, 0.4792890
3.297135, 6.722222, 10.147309, 0, 0.000112, 0.4899760
3.494207, 6.833333, 10.17246, 0, 0.000177, 0.5007130
3.691170, 6.944444, 10.197719, 0, 0.000277, 0.5114950
3.888016, 7.055556, 10.223095, 0, 0.000428, 0.5223200
4.084734, 7.166667, 10.248599, 0, 0.000652, 0.5331820
4.281315, 7.277778, 10.274241, 0, 0.00098, 0.5440770
4.290992, 7.283250, 10.27551, 0, 0.001, 0.5446150
4.477745, 7.388889, 10.300033, 0, 0.001455000, 0.5550040
4.674012, 7.5, 10.325988, 0, 0.00213, 0.5659570
4.870099, 7.611111, 10.352123, 0, 0.003078000, 0.5769360
5.065991, 7.722222, 10.378454, 0, 0.004391000, 0.5879360
5.139437, 7.763912, 10.38839, 0, 0.005, 0.5920690
5.261667, 7.833333, 10.40500, 0, 0.006183000, 0.5989570
5.457105, 7.944444, 10.431784, 0, 0.008595000, 0.6099980
5.549510, 7.997030, 10.44455, 0, 0.01, 0.6152300
5.652280, 8.055556, 10.458831, 1.e-06, 0.01179600, 0.6210570
5.847164, 8.166667, 10.486169, 2.e-06, 0.01598400, 0.6321360
5.996179, 8.251749, 10.50732, 4.e-06, 0.02, 0.6406330
6.041725, 8.277778, 10.513831, 5.e-06, 0.02138800, 0.6432350
6.235924, 8.388889, 10.541854, 1.3e-05, 0.02826100, 0.6543570
6.278641, 8.413360, 10.54808, 1.6e-05, 0.03, 0.6568100
6.429717, 8.5, 10.570283, 3.1e-05, 0.03688200, 0.6655050
6.490555, 8.534934, 10.57931, 4.1e-05, 0.04, 0.6690160
6.623055, 8.611111, 10.599167, 7.3e-05, 0.04754300, 0.6766840
6.662517, 8.633825, 10.60513, 8.7e-05, 0.05, 0.6789740
6.808555, 8.717996, 10.62744, 0.00016, 0.06, 0.6874730
6.815879, 8.722222, 10.628566, 0.000165, 0.06054100, 0.6879010
6.936329, 8.791798, 10.64727, 0.000269, 0.07, 0.6949470
7.008118, 8.833333, 10.658548, 0.000356, 0.07616600, 0.6991630
7.050500, 8.857879, 10.66526, 0.00042, 0.08, 0.7016590
7.154124, 8.917977, 10.68183, 0.000622, 0.09, 0.7077800
7.199693, 8.944444, 10.689196, 0.000737, 0.09468800, 0.7104830
7.249320, 8.973297, 10.69727, 0.000884, 0.1, 0.7134330
7.337620, 9.024712, 10.71180, 0.001215000, 0.11, 0.7187020
7.390505, 9.055556, 10.720606, 0.001463000, 0.1163330, 0.7218710
7.420168, 9.072873, 10.72558, 0.001622000, 0.12, 0.7236530
7.497845, 9.118281, 10.73872, 0.002114000, 0.13, 0.7283370
7.571342, 9.161331, 10.75132, 0.002702000, 0.14, 0.7327920
7.580442, 9.166667, 10.752892, 0.002784000, 0.1412750, 0.7333450
7.641212, 9.202338, 10.76346, 0.003393000, 0.15, 0.7370500
7.707905, 9.241560, 10.77522, 0.004197000, 0.16, 0.7411360
7.769364, 9.277778, 10.786191, 0.005085000, 0.1696120, 0.7449230
7.771792, 9.279210, 10.78663, 0.005123000, 0.17, 0.7450730
7.833186, 9.315465, 10.79775, 0.00618, 0.18, 0.7488780
7.892348, 9.350477, 10.80861, 0.007377000, 0.19, 0.7525660
7.949501, 9.384372, 10.81924, 0.008722000, 0.2, 0.7561500
7.957108, 9.388889, 10.82067, 0.008916000, 0.2013560, 0.7566280
8.004840, 9.417264, 10.82969, 0.01022500, 0.21, 0.7596410
8.058530, 9.449247, 10.83996, 0.01189400, 0.22, 0.7630490
8.110717, 9.480406, 10.85009, 0.01373800, 0.23, 0.7663830
8.143472, 9.5, 10.856528, 0.01501800, 0.2364170, 0.7684870
8.161530, 9.510815, 10.86010, 0.01576600, 0.24, 0.7696500
8.211082, 9.540542, 10.87000, 0.01798500, 0.25, 0.7728580
8.259474, 9.569643, 10.87981, 0.02040400, 0.26, 0.7760110
8.306795, 9.598173, 10.88955, 0.02303000, 0.27, 0.7791170
8.328215, 9.611111, 10.894008, 0.02431000, 0.2745980, 0.7805300
8.353126, 9.626179, 10.89923, 0.02587200, 0.28, 0.7821790
8.398541, 9.653703, 10.90887, 0.02893700, 0.29, 0.7852020
8.443105, 9.680786, 10.91847, 0.03223200, 0.3, 0.7881910
8.486879, 9.707464, 10.92805, 0.03576400, 0.31, 0.7911490
8.511040, 9.722222, 10.933405, 0.03784500, 0.3155940, 0.7927920
8.529915, 9.733769, 10.93762, 0.03953900, 0.32, 0.7940800
8.572265, 9.759732, 10.94720, 0.04356500, 0.33, 0.7969870
8.613974, 9.785382, 10.95679, 0.04784700, 0.34, 0.7998740
8.655085, 9.810744, 10.96640, 0.05239100, 0.35, 0.8027430
8.691586, 9.833333, 10.975081, 0.05670600, 0.3589950, 0.8053110
8.695636, 9.835844, 10.97605, 0.05720200, 0.36, 0.8055970
8.735663, 9.860704, 10.98575, 0.06228500, 0.37, 0.8084400
8.775199, 9.885347, 10.99550, 0.06764600, 0.38, 0.8112720
8.814276, 9.909793, 11.00531, 0.07328800, 0.39, 0.8140970
8.852923, 9.934061, 11.01520, 0.07921400, 0.4, 0.8169180
8.869412, 9.944444, 11.019477, 0.08185100, 0.4042990, 0.8181290
8.891166, 9.958171, 11.02518, 0.08543000, 0.41, 0.8197360
8.929030, 9.982140, 11.03525, 0.09193600, 0.42, 0.8225530
8.966540, 10.005985, 11.04543, 0.09873700, 0.43, 0.8253720
9.003718, 10.029724, 11.05573, 0.1058330, 0.44, 0.8281950
9.040583, 10.053372, 11.06616, 0.1132270, 0.45, 0.8310240
9.043979, 10.055556, 11.067132, 0.1139260, 0.4509250, 0.8312860
9.077157, 10.076945, 11.07673, 0.1209200, 0.46, 0.8338610
9.113456, 10.100458, 11.08746, 0.1289110, 0.47, 0.8367070
9.149500, 10.123927, 11.09835, 0.1372010, 0.48, 0.8395650
9.185304, 10.147366, 11.10943, 0.1457900, 0.49, 0.8424370
9.214635, 10.166667, 11.118699, 0.1530890, 0.4982390, 0.8448150
9.220884, 10.170791, 11.12070, 0.1546750, 0.5, 0.8453250
9.256255, 10.194216, 11.13218, 0.1638560, 0.51, 0.8482300
9.291430, 10.217655, 11.14388, 0.1733300, 0.52, 0.8511540
9.326423, 10.241124, 11.15583, 0.1830940, 0.53, 0.8540990
9.361246, 10.264637, 11.16803, 0.1931430, 0.54, 0.8570680
9.380601, 10.277778, 11.174955, 0.1988720, 0.5455780, 0.8587350
9.395911, 10.28821, 11.18051, 0.2034750, 0.55, 0.8600610
9.430430, 10.311858, 11.19329, 0.2140840, 0.56, 0.8630820
9.464812, 10.335597, 11.20638, 0.2249640, 0.57, 0.8661300
9.499068, 10.359442, 11.21982, 0.2361090, 0.58, 0.8692090
9.533207, 10.383411, 11.23362, 0.2475120, 0.59, 0.8723210
9.540967, 10.388889, 11.236811, 0.2501450, 0.5922770, 0.8730340
9.567238, 10.407521, 11.24780, 0.2591660, 0.6, 0.8754650
9.601170, 10.431789, 11.26241, 0.2710620, 0.61, 0.8786460
9.635010, 10.456235, 11.27746, 0.2831900, 0.62, 0.8818630
9.668765, 10.480878, 11.29299, 0.2955420, 0.63, 0.8851190
9.694703, 10.5, 11.305297, 0.3051990, 0.6377000, 0.8876530
9.702443, 10.505738, 11.30903, 0.3081070, 0.64, 0.8884150
9.736049, 10.530838, 11.32563, 0.3208730, 0.65, 0.8917520
9.769590, 10.55620, 11.34281, 0.3338290, 0.66, 0.8951320
9.803071, 10.58185, 11.36063, 0.3469630, 0.67, 0.8985540
9.836498, 10.607813, 11.37913, 0.3602620, 0.68, 0.9020210
9.840710, 10.611111, 11.381512, 0.3619500, 0.6812610, 0.9024620
9.869875, 10.634118, 11.39836, 0.3737120, 0.69, 0.9055330
9.903207, 10.660796, 11.41838, 0.3872990, 0.7, 0.9090890
9.936498, 10.687879, 11.43926, 0.4010080, 0.71, 0.9126900
9.969753, 10.715403, 11.46105, 0.4148250, 0.72, 0.9163350
9.977900, 10.722222, 11.466545, 0.4182270, 0.7224510, 0.9172350
10.002975, 10.743409, 11.48384, 0.4287350, 0.73, 0.9200220
10.03617, 10.771939, 11.50771, 0.4427220, 0.74, 0.9237510
10.069341, 10.80104, 11.53274, 0.4567710, 0.75, 0.9275180
10.102495, 10.830767, 11.55904, 0.4708670, 0.76, 0.9313200
10.105323, 10.833333, 11.561344, 0.4720710, 0.7608530, 0.9316460
10.135637, 10.861176, 11.58672, 0.4849950, 0.77, 0.9351540
10.168774, 10.892335, 11.61590, 0.4991390, 0.78, 0.9390130
10.201916, 10.924318, 11.64672, 0.5132860, 0.79, 0.9428930
10.222317, 10.944444, 11.666572, 0.5219880, 0.7961530, 0.9452870
10.235074, 10.95721, 11.67934, 0.5274240, 0.8, 0.9467850
10.268263, 10.991105, 11.71395, 0.5415400, 0.81, 0.9506800
10.301499, 11.026116, 11.75073, 0.5556240, 0.82, 0.9545670
10.328631, 11.055556, 11.78248, 0.5670700, 0.8281480, 0.9577200
10.334807, 11.062372, 11.78994, 0.5696680, 0.83, 0.9584350
10.368215, 11.100021, 11.83183, 0.5836670, 0.84, 0.9622680
10.401762, 11.139244, 11.87673, 0.5976170, 0.85, 0.9660510
10.42447, 11.166667, 11.908864, 0.6069910, 0.8567390, 0.9685630
10.435495, 11.180251, 11.92501, 0.6115210, 0.86, 0.9697650
10.469477, 11.223301, 11.97712, 0.6253840, 0.87, 0.9733900
10.50379, 11.268709, 12.03363, 0.6392200, 0.88, 0.9769020
10.510459, 11.277778, 12.045096, 0.6418890, 0.8819300, 0.9775650
10.538537, 11.31687, 12.09520, 0.6530470, 0.89, 0.9802770
10.57386, 11.368284, 12.16271, 0.6668980, 0.9, 0.9834860
10.587526, 11.388889, 12.190252, 0.6721970, 0.9038150, 0.9846610
10.609944, 11.423605, 12.23727, 0.6808150, 0.91, 0.9865000
10.647043, 11.483703, 12.32036, 0.6948630, 0.92, 0.9892890
10.656749, 11.5, 12.343251, 0.6984920, 0.9225610, 0.9899630
10.685514, 11.549784, 12.41405, 0.7091330, 0.93, 0.9918190
10.719233, 11.611111, 12.502989, 0.7213780, 0.9383930, 0.9937180
10.725874, 11.623586, 12.52130, 0.7237600, 0.94, 0.9940570
10.768913, 11.707757, 12.64660, 0.7389480, 0.95, 0.9959710
10.776016, 11.722222, 12.668428, 0.7414130, 0.9515760, 0.9962410
10.815924, 11.806648, 12.79737, 0.7550340, 0.96, 0.9975300
10.828023, 11.833333, 12.838644, 0.7590860, 0.9624000, 0.9978490
10.869256, 11.928222, 12.98719, 0.7726170, 0.97, 0.9987110
10.876046, 11.944444, 13.012843, 0.7748040, 0.9711620, 0.9988230
10.920748, 12.055556, 13.190363, 0.7888970, 0.9781560, 0.9993840
10.933957, 12.089833, 13.24571, 0.7929600, 0.98, 0.9995000
10.96268, 12.166667, 13.370653, 0.8016350, 0.9836590, 0.9996920
11.002292, 12.277778, 13.553264, 0.8132320, 0.9879300, 0.9998530
11.025138, 12.344552, 13.66397, 0.8197250, 0.99, 0.9999070
11.039952, 12.388889, 13.737825, 0.8238590, 0.9911970, 0.9999330
11.075964, 12.5, 13.924036, 0.8336550, 0.9936610, 0.9999700
11.100291, 12.577669, 14.05505, 0.8400700, 0.995, 0.9999840
11.110574, 12.611111, 14.111649, 0.8427320, 0.9954940, 0.9999880
11.143986, 12.722222, 14.300459, 0.8511800, 0.9968380, 0.9999950
11.176368, 12.833333, 14.490299, 0.8590730, 0.9978100, 0.9999980
11.207861, 12.944444, 14.681028, 0.8664720, 0.9985030, 0.9999990
11.238582, 13.055556, 14.872529, 0.8734280, 0.9989900, 1
11.239341, 13.058332, 14.87732, 0.8735960, 0.999, 1
11.268629, 13.166667, 15.064705, 0.8799830, 0.9993270, 1
11.298084, 13.277778, 15.257472, 0.8861730, 0.9995580, 1
11.327018, 13.388889, 15.45076, 0.8920280, 0.9997130, 1
11.35549, 13.5, 15.64451, 0.8975760, 0.9998170, 1
11.383553, 13.611111, 15.838669, 0.9028380, 0.9998840, 1
11.41125, 13.722222, 16.033194, 0.9078340, 0.9999280, 1
11.438619, 13.833333, 16.228047, 0.9125810, 0.9999560, 1
11.465694, 13.944444, 16.423194, 0.9170960, 0.9999730, 1
11.492504, 14.055556, 16.618607, 0.9213910, 0.9999840, 1
11.519074, 14.166667, 16.81426, 0.9254790, 0.9999910, 1
11.545426, 14.277778, 17.01013, 0.9293720, 0.9999940, 1
11.57158, 14.388889, 17.206198, 0.9330780, 0.9999970, 1
11.597553, 14.5, 17.402447, 0.9366090, 0.9999980, 1
11.623361, 14.611111, 17.598862, 0.9399710, 0.9999990, 1
11.649017, 14.722222, 17.795427, 0.9431740, 0.9999990, 1
11.674534, 14.833333, 17.992133, 0.9462240, 1, 1
11.699923, 14.944444, 18.188966, 0.9491290, 1, 1
11.725194, 15.055556, 18.385917, 0.9518950, 1, 1
11.750355, 15.166667, 18.582978, 0.9545280, 1, 1
11.775415, 15.277778, 18.78014, 0.9570340, 1, 1
11.800381, 15.388889, 18.977397, 0.9594190, 1, 1
11.82526, 15.5, 19.17474, 0.9616870, 1, 1
11.850057, 15.611111, 19.372165, 0.9638430, 1, 1
11.874778, 15.722222, 19.569666, 0.9658940, 1, 1
11.899429, 15.833333, 19.767238, 0.9678420, 1, 1
11.924013, 15.944444, 19.964876, 0.9696920, 1, 1
11.948534, 16.055556, 20.162577, 0.9714490, 1, 1
11.972998, 16.166667, 20.360335, 0.9731170, 1, 1
11.997407, 16.277778, 20.558149, 0.9746990, 1, 1
12.021764, 16.388889, 20.756014, 0.9761990, 1, 1
12.046074, 16.5, 20.953926, 0.9776210, 1, 1 ];
%% Matlab GLM fitting
xWT = [5.5,16.5,11,13.8,10.1,14.7,10.4,11.7,9.7,7.3,7.8,8.1,12.2,8.5,11.8,11.7121,11.4083,11.1558,12.4633,12.2761,12.1107,11.9628,11.8291,11.7072,11.5952,11.4917,11.3955,11.3057,11.2214,11.1421]’;
yWT = [0,1,0,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1]’;
n = ones(size(yWT));
[b, dev, stats] = glmfit(xWT, [yWT, n], ‘binomial’, ‘Link’, ‘probit’);
wt = 5 : 0.05 : 17;
[y_fit, dylo, dyhi] = glmval(b, wt, ‘probit’, stats, ‘confidence’, 0.95);
mu_1 = -b(1) /b(2);
gradient = central_diff(y_fit, wt);
go = yWT == 1;
nogo = yWT == 0;
fig1 = figure;
plot(xWT(go), yWT(go), ‘ko’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘g’, ‘MarkerSize’, 8)
hold on
plot(xWT(nogo), yWT(nogo), ‘ks’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘r’, ‘MarkerSize’, 8)
p1 = plot(wt, y_fit, ‘-m’);
p2 = plot(wt, y_fit – dylo, ‘-k’);
plot(wt, y_fit + dyhi, ‘-k’)
p3 = plot(wt, gradient/max(gradient), ‘-‘, ‘Color’, 0.6*[1,1,1]);
xlim([5, 17]);
title(‘Wu and Tian Example’)
xlabel(‘Quantile (q, in units wt)’)
ylabel(‘Probability (p)’)
grid on; grid minor
mfw_pos = get(gcf, ‘Position’); mfw_pos(3) = mfw_pos(3) * 1.6; set(gcf, ‘Position’, mfw_pos); % Make Figure Wide (60% wider)
% These are the results from the R code "Gonogo"
% plot(R_output(:,2), R_output(:,5), ‘–b’)
p4 = plot(R_output(:,2), R_output(:,6), ‘–b’);
plot(R_output(:,2), R_output(:,4), ‘–b’)
p5 = plot(R_output(:,1), R_output(:,5), ‘–r’);
plot(R_output(:,3), R_output(:,5), ‘–r’)
legend([p1 p2 p3 p4 p5],{‘Logistic Fit’,’Confidence Bounds from glmval’,’diff of Logistic Fit’,’Confidence Interval about p’,’Confidence Interval about q’}, …
‘Location’, ‘east’)
mdl = fitglm(xWT, yWT, ‘Distribution’, ‘binomial’, ‘Link’, ‘probit’);
ci = coefCI(mdl, 0.95);
figure(fig1)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,1)), ‘-c’, ‘DisplayName’, ‘from Matlab coefCI’)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,2)), ‘-c’, ‘HandleVisibility’,’off’)I have results from a binary study. The study was analyzied by others using Gonogo (https://www2.isye.gatech.edu/~jeffwu/sensitivitytesting/), the prefered code developed by/for the government for such studies. My code for anlaysis is shown below. I am able to fit a logistic curve using glmfit (magentia line in plot). I then compute the confidence limits about the probablility p using glmval (black lines). These black line match well with the results given by Gonogo (dashed blue lines). However, I am unable to figure out how to compute the confidence limits about the quantile q (the Gonongo results are shown as dashed red lines). I tried using coefCI, with the results in cyan. It is doing something but I don’t understand exactly what, and it does not really match the dashed red lines that I am looking for. Does anyone know how I can compute the confidence limits about the quantile q? Your assistance is appreciated.
%% Output from Gonogo – Main code in next section
R_mu = 10.1707909;
R_sigma = 0.9344091;
R_b0 = -10.88473;
R_b1 = 1.070195;
R_output = [
1.124021, 5.5, 9.875979, 0, 0, 0.3761890
1.321891, 5.611111, 9.900331, 0, 1.e-06, 0.3861210
1.519710, 5.722222, 9.924735, 0, 1.e-06, 0.3961490
1.717474, 5.833333, 9.949193, 0, 2.e-06, 0.4062690
1.915181, 5.944444, 9.973708, 0, 3.e-06, 0.4164760
2.112825, 6.055556, 9.998286, 0, 5.e-06, 0.4267660
2.310404, 6.166667, 10.022929, 0, 9.e-06, 0.4371330
2.507913, 6.277778, 10.047643, 0, 1.5e-05, 0.4475740
2.705346, 6.388889, 10.072432, 0, 2.6e-05, 0.4580840
2.902698, 6.5, 10.097302, 0, 4.3e-05, 0.4686570
3.099963, 6.611111, 10.122259, 0, 7.e-05, 0.4792890
3.297135, 6.722222, 10.147309, 0, 0.000112, 0.4899760
3.494207, 6.833333, 10.17246, 0, 0.000177, 0.5007130
3.691170, 6.944444, 10.197719, 0, 0.000277, 0.5114950
3.888016, 7.055556, 10.223095, 0, 0.000428, 0.5223200
4.084734, 7.166667, 10.248599, 0, 0.000652, 0.5331820
4.281315, 7.277778, 10.274241, 0, 0.00098, 0.5440770
4.290992, 7.283250, 10.27551, 0, 0.001, 0.5446150
4.477745, 7.388889, 10.300033, 0, 0.001455000, 0.5550040
4.674012, 7.5, 10.325988, 0, 0.00213, 0.5659570
4.870099, 7.611111, 10.352123, 0, 0.003078000, 0.5769360
5.065991, 7.722222, 10.378454, 0, 0.004391000, 0.5879360
5.139437, 7.763912, 10.38839, 0, 0.005, 0.5920690
5.261667, 7.833333, 10.40500, 0, 0.006183000, 0.5989570
5.457105, 7.944444, 10.431784, 0, 0.008595000, 0.6099980
5.549510, 7.997030, 10.44455, 0, 0.01, 0.6152300
5.652280, 8.055556, 10.458831, 1.e-06, 0.01179600, 0.6210570
5.847164, 8.166667, 10.486169, 2.e-06, 0.01598400, 0.6321360
5.996179, 8.251749, 10.50732, 4.e-06, 0.02, 0.6406330
6.041725, 8.277778, 10.513831, 5.e-06, 0.02138800, 0.6432350
6.235924, 8.388889, 10.541854, 1.3e-05, 0.02826100, 0.6543570
6.278641, 8.413360, 10.54808, 1.6e-05, 0.03, 0.6568100
6.429717, 8.5, 10.570283, 3.1e-05, 0.03688200, 0.6655050
6.490555, 8.534934, 10.57931, 4.1e-05, 0.04, 0.6690160
6.623055, 8.611111, 10.599167, 7.3e-05, 0.04754300, 0.6766840
6.662517, 8.633825, 10.60513, 8.7e-05, 0.05, 0.6789740
6.808555, 8.717996, 10.62744, 0.00016, 0.06, 0.6874730
6.815879, 8.722222, 10.628566, 0.000165, 0.06054100, 0.6879010
6.936329, 8.791798, 10.64727, 0.000269, 0.07, 0.6949470
7.008118, 8.833333, 10.658548, 0.000356, 0.07616600, 0.6991630
7.050500, 8.857879, 10.66526, 0.00042, 0.08, 0.7016590
7.154124, 8.917977, 10.68183, 0.000622, 0.09, 0.7077800
7.199693, 8.944444, 10.689196, 0.000737, 0.09468800, 0.7104830
7.249320, 8.973297, 10.69727, 0.000884, 0.1, 0.7134330
7.337620, 9.024712, 10.71180, 0.001215000, 0.11, 0.7187020
7.390505, 9.055556, 10.720606, 0.001463000, 0.1163330, 0.7218710
7.420168, 9.072873, 10.72558, 0.001622000, 0.12, 0.7236530
7.497845, 9.118281, 10.73872, 0.002114000, 0.13, 0.7283370
7.571342, 9.161331, 10.75132, 0.002702000, 0.14, 0.7327920
7.580442, 9.166667, 10.752892, 0.002784000, 0.1412750, 0.7333450
7.641212, 9.202338, 10.76346, 0.003393000, 0.15, 0.7370500
7.707905, 9.241560, 10.77522, 0.004197000, 0.16, 0.7411360
7.769364, 9.277778, 10.786191, 0.005085000, 0.1696120, 0.7449230
7.771792, 9.279210, 10.78663, 0.005123000, 0.17, 0.7450730
7.833186, 9.315465, 10.79775, 0.00618, 0.18, 0.7488780
7.892348, 9.350477, 10.80861, 0.007377000, 0.19, 0.7525660
7.949501, 9.384372, 10.81924, 0.008722000, 0.2, 0.7561500
7.957108, 9.388889, 10.82067, 0.008916000, 0.2013560, 0.7566280
8.004840, 9.417264, 10.82969, 0.01022500, 0.21, 0.7596410
8.058530, 9.449247, 10.83996, 0.01189400, 0.22, 0.7630490
8.110717, 9.480406, 10.85009, 0.01373800, 0.23, 0.7663830
8.143472, 9.5, 10.856528, 0.01501800, 0.2364170, 0.7684870
8.161530, 9.510815, 10.86010, 0.01576600, 0.24, 0.7696500
8.211082, 9.540542, 10.87000, 0.01798500, 0.25, 0.7728580
8.259474, 9.569643, 10.87981, 0.02040400, 0.26, 0.7760110
8.306795, 9.598173, 10.88955, 0.02303000, 0.27, 0.7791170
8.328215, 9.611111, 10.894008, 0.02431000, 0.2745980, 0.7805300
8.353126, 9.626179, 10.89923, 0.02587200, 0.28, 0.7821790
8.398541, 9.653703, 10.90887, 0.02893700, 0.29, 0.7852020
8.443105, 9.680786, 10.91847, 0.03223200, 0.3, 0.7881910
8.486879, 9.707464, 10.92805, 0.03576400, 0.31, 0.7911490
8.511040, 9.722222, 10.933405, 0.03784500, 0.3155940, 0.7927920
8.529915, 9.733769, 10.93762, 0.03953900, 0.32, 0.7940800
8.572265, 9.759732, 10.94720, 0.04356500, 0.33, 0.7969870
8.613974, 9.785382, 10.95679, 0.04784700, 0.34, 0.7998740
8.655085, 9.810744, 10.96640, 0.05239100, 0.35, 0.8027430
8.691586, 9.833333, 10.975081, 0.05670600, 0.3589950, 0.8053110
8.695636, 9.835844, 10.97605, 0.05720200, 0.36, 0.8055970
8.735663, 9.860704, 10.98575, 0.06228500, 0.37, 0.8084400
8.775199, 9.885347, 10.99550, 0.06764600, 0.38, 0.8112720
8.814276, 9.909793, 11.00531, 0.07328800, 0.39, 0.8140970
8.852923, 9.934061, 11.01520, 0.07921400, 0.4, 0.8169180
8.869412, 9.944444, 11.019477, 0.08185100, 0.4042990, 0.8181290
8.891166, 9.958171, 11.02518, 0.08543000, 0.41, 0.8197360
8.929030, 9.982140, 11.03525, 0.09193600, 0.42, 0.8225530
8.966540, 10.005985, 11.04543, 0.09873700, 0.43, 0.8253720
9.003718, 10.029724, 11.05573, 0.1058330, 0.44, 0.8281950
9.040583, 10.053372, 11.06616, 0.1132270, 0.45, 0.8310240
9.043979, 10.055556, 11.067132, 0.1139260, 0.4509250, 0.8312860
9.077157, 10.076945, 11.07673, 0.1209200, 0.46, 0.8338610
9.113456, 10.100458, 11.08746, 0.1289110, 0.47, 0.8367070
9.149500, 10.123927, 11.09835, 0.1372010, 0.48, 0.8395650
9.185304, 10.147366, 11.10943, 0.1457900, 0.49, 0.8424370
9.214635, 10.166667, 11.118699, 0.1530890, 0.4982390, 0.8448150
9.220884, 10.170791, 11.12070, 0.1546750, 0.5, 0.8453250
9.256255, 10.194216, 11.13218, 0.1638560, 0.51, 0.8482300
9.291430, 10.217655, 11.14388, 0.1733300, 0.52, 0.8511540
9.326423, 10.241124, 11.15583, 0.1830940, 0.53, 0.8540990
9.361246, 10.264637, 11.16803, 0.1931430, 0.54, 0.8570680
9.380601, 10.277778, 11.174955, 0.1988720, 0.5455780, 0.8587350
9.395911, 10.28821, 11.18051, 0.2034750, 0.55, 0.8600610
9.430430, 10.311858, 11.19329, 0.2140840, 0.56, 0.8630820
9.464812, 10.335597, 11.20638, 0.2249640, 0.57, 0.8661300
9.499068, 10.359442, 11.21982, 0.2361090, 0.58, 0.8692090
9.533207, 10.383411, 11.23362, 0.2475120, 0.59, 0.8723210
9.540967, 10.388889, 11.236811, 0.2501450, 0.5922770, 0.8730340
9.567238, 10.407521, 11.24780, 0.2591660, 0.6, 0.8754650
9.601170, 10.431789, 11.26241, 0.2710620, 0.61, 0.8786460
9.635010, 10.456235, 11.27746, 0.2831900, 0.62, 0.8818630
9.668765, 10.480878, 11.29299, 0.2955420, 0.63, 0.8851190
9.694703, 10.5, 11.305297, 0.3051990, 0.6377000, 0.8876530
9.702443, 10.505738, 11.30903, 0.3081070, 0.64, 0.8884150
9.736049, 10.530838, 11.32563, 0.3208730, 0.65, 0.8917520
9.769590, 10.55620, 11.34281, 0.3338290, 0.66, 0.8951320
9.803071, 10.58185, 11.36063, 0.3469630, 0.67, 0.8985540
9.836498, 10.607813, 11.37913, 0.3602620, 0.68, 0.9020210
9.840710, 10.611111, 11.381512, 0.3619500, 0.6812610, 0.9024620
9.869875, 10.634118, 11.39836, 0.3737120, 0.69, 0.9055330
9.903207, 10.660796, 11.41838, 0.3872990, 0.7, 0.9090890
9.936498, 10.687879, 11.43926, 0.4010080, 0.71, 0.9126900
9.969753, 10.715403, 11.46105, 0.4148250, 0.72, 0.9163350
9.977900, 10.722222, 11.466545, 0.4182270, 0.7224510, 0.9172350
10.002975, 10.743409, 11.48384, 0.4287350, 0.73, 0.9200220
10.03617, 10.771939, 11.50771, 0.4427220, 0.74, 0.9237510
10.069341, 10.80104, 11.53274, 0.4567710, 0.75, 0.9275180
10.102495, 10.830767, 11.55904, 0.4708670, 0.76, 0.9313200
10.105323, 10.833333, 11.561344, 0.4720710, 0.7608530, 0.9316460
10.135637, 10.861176, 11.58672, 0.4849950, 0.77, 0.9351540
10.168774, 10.892335, 11.61590, 0.4991390, 0.78, 0.9390130
10.201916, 10.924318, 11.64672, 0.5132860, 0.79, 0.9428930
10.222317, 10.944444, 11.666572, 0.5219880, 0.7961530, 0.9452870
10.235074, 10.95721, 11.67934, 0.5274240, 0.8, 0.9467850
10.268263, 10.991105, 11.71395, 0.5415400, 0.81, 0.9506800
10.301499, 11.026116, 11.75073, 0.5556240, 0.82, 0.9545670
10.328631, 11.055556, 11.78248, 0.5670700, 0.8281480, 0.9577200
10.334807, 11.062372, 11.78994, 0.5696680, 0.83, 0.9584350
10.368215, 11.100021, 11.83183, 0.5836670, 0.84, 0.9622680
10.401762, 11.139244, 11.87673, 0.5976170, 0.85, 0.9660510
10.42447, 11.166667, 11.908864, 0.6069910, 0.8567390, 0.9685630
10.435495, 11.180251, 11.92501, 0.6115210, 0.86, 0.9697650
10.469477, 11.223301, 11.97712, 0.6253840, 0.87, 0.9733900
10.50379, 11.268709, 12.03363, 0.6392200, 0.88, 0.9769020
10.510459, 11.277778, 12.045096, 0.6418890, 0.8819300, 0.9775650
10.538537, 11.31687, 12.09520, 0.6530470, 0.89, 0.9802770
10.57386, 11.368284, 12.16271, 0.6668980, 0.9, 0.9834860
10.587526, 11.388889, 12.190252, 0.6721970, 0.9038150, 0.9846610
10.609944, 11.423605, 12.23727, 0.6808150, 0.91, 0.9865000
10.647043, 11.483703, 12.32036, 0.6948630, 0.92, 0.9892890
10.656749, 11.5, 12.343251, 0.6984920, 0.9225610, 0.9899630
10.685514, 11.549784, 12.41405, 0.7091330, 0.93, 0.9918190
10.719233, 11.611111, 12.502989, 0.7213780, 0.9383930, 0.9937180
10.725874, 11.623586, 12.52130, 0.7237600, 0.94, 0.9940570
10.768913, 11.707757, 12.64660, 0.7389480, 0.95, 0.9959710
10.776016, 11.722222, 12.668428, 0.7414130, 0.9515760, 0.9962410
10.815924, 11.806648, 12.79737, 0.7550340, 0.96, 0.9975300
10.828023, 11.833333, 12.838644, 0.7590860, 0.9624000, 0.9978490
10.869256, 11.928222, 12.98719, 0.7726170, 0.97, 0.9987110
10.876046, 11.944444, 13.012843, 0.7748040, 0.9711620, 0.9988230
10.920748, 12.055556, 13.190363, 0.7888970, 0.9781560, 0.9993840
10.933957, 12.089833, 13.24571, 0.7929600, 0.98, 0.9995000
10.96268, 12.166667, 13.370653, 0.8016350, 0.9836590, 0.9996920
11.002292, 12.277778, 13.553264, 0.8132320, 0.9879300, 0.9998530
11.025138, 12.344552, 13.66397, 0.8197250, 0.99, 0.9999070
11.039952, 12.388889, 13.737825, 0.8238590, 0.9911970, 0.9999330
11.075964, 12.5, 13.924036, 0.8336550, 0.9936610, 0.9999700
11.100291, 12.577669, 14.05505, 0.8400700, 0.995, 0.9999840
11.110574, 12.611111, 14.111649, 0.8427320, 0.9954940, 0.9999880
11.143986, 12.722222, 14.300459, 0.8511800, 0.9968380, 0.9999950
11.176368, 12.833333, 14.490299, 0.8590730, 0.9978100, 0.9999980
11.207861, 12.944444, 14.681028, 0.8664720, 0.9985030, 0.9999990
11.238582, 13.055556, 14.872529, 0.8734280, 0.9989900, 1
11.239341, 13.058332, 14.87732, 0.8735960, 0.999, 1
11.268629, 13.166667, 15.064705, 0.8799830, 0.9993270, 1
11.298084, 13.277778, 15.257472, 0.8861730, 0.9995580, 1
11.327018, 13.388889, 15.45076, 0.8920280, 0.9997130, 1
11.35549, 13.5, 15.64451, 0.8975760, 0.9998170, 1
11.383553, 13.611111, 15.838669, 0.9028380, 0.9998840, 1
11.41125, 13.722222, 16.033194, 0.9078340, 0.9999280, 1
11.438619, 13.833333, 16.228047, 0.9125810, 0.9999560, 1
11.465694, 13.944444, 16.423194, 0.9170960, 0.9999730, 1
11.492504, 14.055556, 16.618607, 0.9213910, 0.9999840, 1
11.519074, 14.166667, 16.81426, 0.9254790, 0.9999910, 1
11.545426, 14.277778, 17.01013, 0.9293720, 0.9999940, 1
11.57158, 14.388889, 17.206198, 0.9330780, 0.9999970, 1
11.597553, 14.5, 17.402447, 0.9366090, 0.9999980, 1
11.623361, 14.611111, 17.598862, 0.9399710, 0.9999990, 1
11.649017, 14.722222, 17.795427, 0.9431740, 0.9999990, 1
11.674534, 14.833333, 17.992133, 0.9462240, 1, 1
11.699923, 14.944444, 18.188966, 0.9491290, 1, 1
11.725194, 15.055556, 18.385917, 0.9518950, 1, 1
11.750355, 15.166667, 18.582978, 0.9545280, 1, 1
11.775415, 15.277778, 18.78014, 0.9570340, 1, 1
11.800381, 15.388889, 18.977397, 0.9594190, 1, 1
11.82526, 15.5, 19.17474, 0.9616870, 1, 1
11.850057, 15.611111, 19.372165, 0.9638430, 1, 1
11.874778, 15.722222, 19.569666, 0.9658940, 1, 1
11.899429, 15.833333, 19.767238, 0.9678420, 1, 1
11.924013, 15.944444, 19.964876, 0.9696920, 1, 1
11.948534, 16.055556, 20.162577, 0.9714490, 1, 1
11.972998, 16.166667, 20.360335, 0.9731170, 1, 1
11.997407, 16.277778, 20.558149, 0.9746990, 1, 1
12.021764, 16.388889, 20.756014, 0.9761990, 1, 1
12.046074, 16.5, 20.953926, 0.9776210, 1, 1 ];
%% Matlab GLM fitting
xWT = [5.5,16.5,11,13.8,10.1,14.7,10.4,11.7,9.7,7.3,7.8,8.1,12.2,8.5,11.8,11.7121,11.4083,11.1558,12.4633,12.2761,12.1107,11.9628,11.8291,11.7072,11.5952,11.4917,11.3955,11.3057,11.2214,11.1421]’;
yWT = [0,1,0,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1]’;
n = ones(size(yWT));
[b, dev, stats] = glmfit(xWT, [yWT, n], ‘binomial’, ‘Link’, ‘probit’);
wt = 5 : 0.05 : 17;
[y_fit, dylo, dyhi] = glmval(b, wt, ‘probit’, stats, ‘confidence’, 0.95);
mu_1 = -b(1) /b(2);
gradient = central_diff(y_fit, wt);
go = yWT == 1;
nogo = yWT == 0;
fig1 = figure;
plot(xWT(go), yWT(go), ‘ko’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘g’, ‘MarkerSize’, 8)
hold on
plot(xWT(nogo), yWT(nogo), ‘ks’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘r’, ‘MarkerSize’, 8)
p1 = plot(wt, y_fit, ‘-m’);
p2 = plot(wt, y_fit – dylo, ‘-k’);
plot(wt, y_fit + dyhi, ‘-k’)
p3 = plot(wt, gradient/max(gradient), ‘-‘, ‘Color’, 0.6*[1,1,1]);
xlim([5, 17]);
title(‘Wu and Tian Example’)
xlabel(‘Quantile (q, in units wt)’)
ylabel(‘Probability (p)’)
grid on; grid minor
mfw_pos = get(gcf, ‘Position’); mfw_pos(3) = mfw_pos(3) * 1.6; set(gcf, ‘Position’, mfw_pos); % Make Figure Wide (60% wider)
% These are the results from the R code "Gonogo"
% plot(R_output(:,2), R_output(:,5), ‘–b’)
p4 = plot(R_output(:,2), R_output(:,6), ‘–b’);
plot(R_output(:,2), R_output(:,4), ‘–b’)
p5 = plot(R_output(:,1), R_output(:,5), ‘–r’);
plot(R_output(:,3), R_output(:,5), ‘–r’)
legend([p1 p2 p3 p4 p5],{‘Logistic Fit’,’Confidence Bounds from glmval’,’diff of Logistic Fit’,’Confidence Interval about p’,’Confidence Interval about q’}, …
‘Location’, ‘east’)
mdl = fitglm(xWT, yWT, ‘Distribution’, ‘binomial’, ‘Link’, ‘probit’);
ci = coefCI(mdl, 0.95);
figure(fig1)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,1)), ‘-c’, ‘DisplayName’, ‘from Matlab coefCI’)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,2)), ‘-c’, ‘HandleVisibility’,’off’) I have results from a binary study. The study was analyzied by others using Gonogo (https://www2.isye.gatech.edu/~jeffwu/sensitivitytesting/), the prefered code developed by/for the government for such studies. My code for anlaysis is shown below. I am able to fit a logistic curve using glmfit (magentia line in plot). I then compute the confidence limits about the probablility p using glmval (black lines). These black line match well with the results given by Gonogo (dashed blue lines). However, I am unable to figure out how to compute the confidence limits about the quantile q (the Gonongo results are shown as dashed red lines). I tried using coefCI, with the results in cyan. It is doing something but I don’t understand exactly what, and it does not really match the dashed red lines that I am looking for. Does anyone know how I can compute the confidence limits about the quantile q? Your assistance is appreciated.
%% Output from Gonogo – Main code in next section
R_mu = 10.1707909;
R_sigma = 0.9344091;
R_b0 = -10.88473;
R_b1 = 1.070195;
R_output = [
1.124021, 5.5, 9.875979, 0, 0, 0.3761890
1.321891, 5.611111, 9.900331, 0, 1.e-06, 0.3861210
1.519710, 5.722222, 9.924735, 0, 1.e-06, 0.3961490
1.717474, 5.833333, 9.949193, 0, 2.e-06, 0.4062690
1.915181, 5.944444, 9.973708, 0, 3.e-06, 0.4164760
2.112825, 6.055556, 9.998286, 0, 5.e-06, 0.4267660
2.310404, 6.166667, 10.022929, 0, 9.e-06, 0.4371330
2.507913, 6.277778, 10.047643, 0, 1.5e-05, 0.4475740
2.705346, 6.388889, 10.072432, 0, 2.6e-05, 0.4580840
2.902698, 6.5, 10.097302, 0, 4.3e-05, 0.4686570
3.099963, 6.611111, 10.122259, 0, 7.e-05, 0.4792890
3.297135, 6.722222, 10.147309, 0, 0.000112, 0.4899760
3.494207, 6.833333, 10.17246, 0, 0.000177, 0.5007130
3.691170, 6.944444, 10.197719, 0, 0.000277, 0.5114950
3.888016, 7.055556, 10.223095, 0, 0.000428, 0.5223200
4.084734, 7.166667, 10.248599, 0, 0.000652, 0.5331820
4.281315, 7.277778, 10.274241, 0, 0.00098, 0.5440770
4.290992, 7.283250, 10.27551, 0, 0.001, 0.5446150
4.477745, 7.388889, 10.300033, 0, 0.001455000, 0.5550040
4.674012, 7.5, 10.325988, 0, 0.00213, 0.5659570
4.870099, 7.611111, 10.352123, 0, 0.003078000, 0.5769360
5.065991, 7.722222, 10.378454, 0, 0.004391000, 0.5879360
5.139437, 7.763912, 10.38839, 0, 0.005, 0.5920690
5.261667, 7.833333, 10.40500, 0, 0.006183000, 0.5989570
5.457105, 7.944444, 10.431784, 0, 0.008595000, 0.6099980
5.549510, 7.997030, 10.44455, 0, 0.01, 0.6152300
5.652280, 8.055556, 10.458831, 1.e-06, 0.01179600, 0.6210570
5.847164, 8.166667, 10.486169, 2.e-06, 0.01598400, 0.6321360
5.996179, 8.251749, 10.50732, 4.e-06, 0.02, 0.6406330
6.041725, 8.277778, 10.513831, 5.e-06, 0.02138800, 0.6432350
6.235924, 8.388889, 10.541854, 1.3e-05, 0.02826100, 0.6543570
6.278641, 8.413360, 10.54808, 1.6e-05, 0.03, 0.6568100
6.429717, 8.5, 10.570283, 3.1e-05, 0.03688200, 0.6655050
6.490555, 8.534934, 10.57931, 4.1e-05, 0.04, 0.6690160
6.623055, 8.611111, 10.599167, 7.3e-05, 0.04754300, 0.6766840
6.662517, 8.633825, 10.60513, 8.7e-05, 0.05, 0.6789740
6.808555, 8.717996, 10.62744, 0.00016, 0.06, 0.6874730
6.815879, 8.722222, 10.628566, 0.000165, 0.06054100, 0.6879010
6.936329, 8.791798, 10.64727, 0.000269, 0.07, 0.6949470
7.008118, 8.833333, 10.658548, 0.000356, 0.07616600, 0.6991630
7.050500, 8.857879, 10.66526, 0.00042, 0.08, 0.7016590
7.154124, 8.917977, 10.68183, 0.000622, 0.09, 0.7077800
7.199693, 8.944444, 10.689196, 0.000737, 0.09468800, 0.7104830
7.249320, 8.973297, 10.69727, 0.000884, 0.1, 0.7134330
7.337620, 9.024712, 10.71180, 0.001215000, 0.11, 0.7187020
7.390505, 9.055556, 10.720606, 0.001463000, 0.1163330, 0.7218710
7.420168, 9.072873, 10.72558, 0.001622000, 0.12, 0.7236530
7.497845, 9.118281, 10.73872, 0.002114000, 0.13, 0.7283370
7.571342, 9.161331, 10.75132, 0.002702000, 0.14, 0.7327920
7.580442, 9.166667, 10.752892, 0.002784000, 0.1412750, 0.7333450
7.641212, 9.202338, 10.76346, 0.003393000, 0.15, 0.7370500
7.707905, 9.241560, 10.77522, 0.004197000, 0.16, 0.7411360
7.769364, 9.277778, 10.786191, 0.005085000, 0.1696120, 0.7449230
7.771792, 9.279210, 10.78663, 0.005123000, 0.17, 0.7450730
7.833186, 9.315465, 10.79775, 0.00618, 0.18, 0.7488780
7.892348, 9.350477, 10.80861, 0.007377000, 0.19, 0.7525660
7.949501, 9.384372, 10.81924, 0.008722000, 0.2, 0.7561500
7.957108, 9.388889, 10.82067, 0.008916000, 0.2013560, 0.7566280
8.004840, 9.417264, 10.82969, 0.01022500, 0.21, 0.7596410
8.058530, 9.449247, 10.83996, 0.01189400, 0.22, 0.7630490
8.110717, 9.480406, 10.85009, 0.01373800, 0.23, 0.7663830
8.143472, 9.5, 10.856528, 0.01501800, 0.2364170, 0.7684870
8.161530, 9.510815, 10.86010, 0.01576600, 0.24, 0.7696500
8.211082, 9.540542, 10.87000, 0.01798500, 0.25, 0.7728580
8.259474, 9.569643, 10.87981, 0.02040400, 0.26, 0.7760110
8.306795, 9.598173, 10.88955, 0.02303000, 0.27, 0.7791170
8.328215, 9.611111, 10.894008, 0.02431000, 0.2745980, 0.7805300
8.353126, 9.626179, 10.89923, 0.02587200, 0.28, 0.7821790
8.398541, 9.653703, 10.90887, 0.02893700, 0.29, 0.7852020
8.443105, 9.680786, 10.91847, 0.03223200, 0.3, 0.7881910
8.486879, 9.707464, 10.92805, 0.03576400, 0.31, 0.7911490
8.511040, 9.722222, 10.933405, 0.03784500, 0.3155940, 0.7927920
8.529915, 9.733769, 10.93762, 0.03953900, 0.32, 0.7940800
8.572265, 9.759732, 10.94720, 0.04356500, 0.33, 0.7969870
8.613974, 9.785382, 10.95679, 0.04784700, 0.34, 0.7998740
8.655085, 9.810744, 10.96640, 0.05239100, 0.35, 0.8027430
8.691586, 9.833333, 10.975081, 0.05670600, 0.3589950, 0.8053110
8.695636, 9.835844, 10.97605, 0.05720200, 0.36, 0.8055970
8.735663, 9.860704, 10.98575, 0.06228500, 0.37, 0.8084400
8.775199, 9.885347, 10.99550, 0.06764600, 0.38, 0.8112720
8.814276, 9.909793, 11.00531, 0.07328800, 0.39, 0.8140970
8.852923, 9.934061, 11.01520, 0.07921400, 0.4, 0.8169180
8.869412, 9.944444, 11.019477, 0.08185100, 0.4042990, 0.8181290
8.891166, 9.958171, 11.02518, 0.08543000, 0.41, 0.8197360
8.929030, 9.982140, 11.03525, 0.09193600, 0.42, 0.8225530
8.966540, 10.005985, 11.04543, 0.09873700, 0.43, 0.8253720
9.003718, 10.029724, 11.05573, 0.1058330, 0.44, 0.8281950
9.040583, 10.053372, 11.06616, 0.1132270, 0.45, 0.8310240
9.043979, 10.055556, 11.067132, 0.1139260, 0.4509250, 0.8312860
9.077157, 10.076945, 11.07673, 0.1209200, 0.46, 0.8338610
9.113456, 10.100458, 11.08746, 0.1289110, 0.47, 0.8367070
9.149500, 10.123927, 11.09835, 0.1372010, 0.48, 0.8395650
9.185304, 10.147366, 11.10943, 0.1457900, 0.49, 0.8424370
9.214635, 10.166667, 11.118699, 0.1530890, 0.4982390, 0.8448150
9.220884, 10.170791, 11.12070, 0.1546750, 0.5, 0.8453250
9.256255, 10.194216, 11.13218, 0.1638560, 0.51, 0.8482300
9.291430, 10.217655, 11.14388, 0.1733300, 0.52, 0.8511540
9.326423, 10.241124, 11.15583, 0.1830940, 0.53, 0.8540990
9.361246, 10.264637, 11.16803, 0.1931430, 0.54, 0.8570680
9.380601, 10.277778, 11.174955, 0.1988720, 0.5455780, 0.8587350
9.395911, 10.28821, 11.18051, 0.2034750, 0.55, 0.8600610
9.430430, 10.311858, 11.19329, 0.2140840, 0.56, 0.8630820
9.464812, 10.335597, 11.20638, 0.2249640, 0.57, 0.8661300
9.499068, 10.359442, 11.21982, 0.2361090, 0.58, 0.8692090
9.533207, 10.383411, 11.23362, 0.2475120, 0.59, 0.8723210
9.540967, 10.388889, 11.236811, 0.2501450, 0.5922770, 0.8730340
9.567238, 10.407521, 11.24780, 0.2591660, 0.6, 0.8754650
9.601170, 10.431789, 11.26241, 0.2710620, 0.61, 0.8786460
9.635010, 10.456235, 11.27746, 0.2831900, 0.62, 0.8818630
9.668765, 10.480878, 11.29299, 0.2955420, 0.63, 0.8851190
9.694703, 10.5, 11.305297, 0.3051990, 0.6377000, 0.8876530
9.702443, 10.505738, 11.30903, 0.3081070, 0.64, 0.8884150
9.736049, 10.530838, 11.32563, 0.3208730, 0.65, 0.8917520
9.769590, 10.55620, 11.34281, 0.3338290, 0.66, 0.8951320
9.803071, 10.58185, 11.36063, 0.3469630, 0.67, 0.8985540
9.836498, 10.607813, 11.37913, 0.3602620, 0.68, 0.9020210
9.840710, 10.611111, 11.381512, 0.3619500, 0.6812610, 0.9024620
9.869875, 10.634118, 11.39836, 0.3737120, 0.69, 0.9055330
9.903207, 10.660796, 11.41838, 0.3872990, 0.7, 0.9090890
9.936498, 10.687879, 11.43926, 0.4010080, 0.71, 0.9126900
9.969753, 10.715403, 11.46105, 0.4148250, 0.72, 0.9163350
9.977900, 10.722222, 11.466545, 0.4182270, 0.7224510, 0.9172350
10.002975, 10.743409, 11.48384, 0.4287350, 0.73, 0.9200220
10.03617, 10.771939, 11.50771, 0.4427220, 0.74, 0.9237510
10.069341, 10.80104, 11.53274, 0.4567710, 0.75, 0.9275180
10.102495, 10.830767, 11.55904, 0.4708670, 0.76, 0.9313200
10.105323, 10.833333, 11.561344, 0.4720710, 0.7608530, 0.9316460
10.135637, 10.861176, 11.58672, 0.4849950, 0.77, 0.9351540
10.168774, 10.892335, 11.61590, 0.4991390, 0.78, 0.9390130
10.201916, 10.924318, 11.64672, 0.5132860, 0.79, 0.9428930
10.222317, 10.944444, 11.666572, 0.5219880, 0.7961530, 0.9452870
10.235074, 10.95721, 11.67934, 0.5274240, 0.8, 0.9467850
10.268263, 10.991105, 11.71395, 0.5415400, 0.81, 0.9506800
10.301499, 11.026116, 11.75073, 0.5556240, 0.82, 0.9545670
10.328631, 11.055556, 11.78248, 0.5670700, 0.8281480, 0.9577200
10.334807, 11.062372, 11.78994, 0.5696680, 0.83, 0.9584350
10.368215, 11.100021, 11.83183, 0.5836670, 0.84, 0.9622680
10.401762, 11.139244, 11.87673, 0.5976170, 0.85, 0.9660510
10.42447, 11.166667, 11.908864, 0.6069910, 0.8567390, 0.9685630
10.435495, 11.180251, 11.92501, 0.6115210, 0.86, 0.9697650
10.469477, 11.223301, 11.97712, 0.6253840, 0.87, 0.9733900
10.50379, 11.268709, 12.03363, 0.6392200, 0.88, 0.9769020
10.510459, 11.277778, 12.045096, 0.6418890, 0.8819300, 0.9775650
10.538537, 11.31687, 12.09520, 0.6530470, 0.89, 0.9802770
10.57386, 11.368284, 12.16271, 0.6668980, 0.9, 0.9834860
10.587526, 11.388889, 12.190252, 0.6721970, 0.9038150, 0.9846610
10.609944, 11.423605, 12.23727, 0.6808150, 0.91, 0.9865000
10.647043, 11.483703, 12.32036, 0.6948630, 0.92, 0.9892890
10.656749, 11.5, 12.343251, 0.6984920, 0.9225610, 0.9899630
10.685514, 11.549784, 12.41405, 0.7091330, 0.93, 0.9918190
10.719233, 11.611111, 12.502989, 0.7213780, 0.9383930, 0.9937180
10.725874, 11.623586, 12.52130, 0.7237600, 0.94, 0.9940570
10.768913, 11.707757, 12.64660, 0.7389480, 0.95, 0.9959710
10.776016, 11.722222, 12.668428, 0.7414130, 0.9515760, 0.9962410
10.815924, 11.806648, 12.79737, 0.7550340, 0.96, 0.9975300
10.828023, 11.833333, 12.838644, 0.7590860, 0.9624000, 0.9978490
10.869256, 11.928222, 12.98719, 0.7726170, 0.97, 0.9987110
10.876046, 11.944444, 13.012843, 0.7748040, 0.9711620, 0.9988230
10.920748, 12.055556, 13.190363, 0.7888970, 0.9781560, 0.9993840
10.933957, 12.089833, 13.24571, 0.7929600, 0.98, 0.9995000
10.96268, 12.166667, 13.370653, 0.8016350, 0.9836590, 0.9996920
11.002292, 12.277778, 13.553264, 0.8132320, 0.9879300, 0.9998530
11.025138, 12.344552, 13.66397, 0.8197250, 0.99, 0.9999070
11.039952, 12.388889, 13.737825, 0.8238590, 0.9911970, 0.9999330
11.075964, 12.5, 13.924036, 0.8336550, 0.9936610, 0.9999700
11.100291, 12.577669, 14.05505, 0.8400700, 0.995, 0.9999840
11.110574, 12.611111, 14.111649, 0.8427320, 0.9954940, 0.9999880
11.143986, 12.722222, 14.300459, 0.8511800, 0.9968380, 0.9999950
11.176368, 12.833333, 14.490299, 0.8590730, 0.9978100, 0.9999980
11.207861, 12.944444, 14.681028, 0.8664720, 0.9985030, 0.9999990
11.238582, 13.055556, 14.872529, 0.8734280, 0.9989900, 1
11.239341, 13.058332, 14.87732, 0.8735960, 0.999, 1
11.268629, 13.166667, 15.064705, 0.8799830, 0.9993270, 1
11.298084, 13.277778, 15.257472, 0.8861730, 0.9995580, 1
11.327018, 13.388889, 15.45076, 0.8920280, 0.9997130, 1
11.35549, 13.5, 15.64451, 0.8975760, 0.9998170, 1
11.383553, 13.611111, 15.838669, 0.9028380, 0.9998840, 1
11.41125, 13.722222, 16.033194, 0.9078340, 0.9999280, 1
11.438619, 13.833333, 16.228047, 0.9125810, 0.9999560, 1
11.465694, 13.944444, 16.423194, 0.9170960, 0.9999730, 1
11.492504, 14.055556, 16.618607, 0.9213910, 0.9999840, 1
11.519074, 14.166667, 16.81426, 0.9254790, 0.9999910, 1
11.545426, 14.277778, 17.01013, 0.9293720, 0.9999940, 1
11.57158, 14.388889, 17.206198, 0.9330780, 0.9999970, 1
11.597553, 14.5, 17.402447, 0.9366090, 0.9999980, 1
11.623361, 14.611111, 17.598862, 0.9399710, 0.9999990, 1
11.649017, 14.722222, 17.795427, 0.9431740, 0.9999990, 1
11.674534, 14.833333, 17.992133, 0.9462240, 1, 1
11.699923, 14.944444, 18.188966, 0.9491290, 1, 1
11.725194, 15.055556, 18.385917, 0.9518950, 1, 1
11.750355, 15.166667, 18.582978, 0.9545280, 1, 1
11.775415, 15.277778, 18.78014, 0.9570340, 1, 1
11.800381, 15.388889, 18.977397, 0.9594190, 1, 1
11.82526, 15.5, 19.17474, 0.9616870, 1, 1
11.850057, 15.611111, 19.372165, 0.9638430, 1, 1
11.874778, 15.722222, 19.569666, 0.9658940, 1, 1
11.899429, 15.833333, 19.767238, 0.9678420, 1, 1
11.924013, 15.944444, 19.964876, 0.9696920, 1, 1
11.948534, 16.055556, 20.162577, 0.9714490, 1, 1
11.972998, 16.166667, 20.360335, 0.9731170, 1, 1
11.997407, 16.277778, 20.558149, 0.9746990, 1, 1
12.021764, 16.388889, 20.756014, 0.9761990, 1, 1
12.046074, 16.5, 20.953926, 0.9776210, 1, 1 ];
%% Matlab GLM fitting
xWT = [5.5,16.5,11,13.8,10.1,14.7,10.4,11.7,9.7,7.3,7.8,8.1,12.2,8.5,11.8,11.7121,11.4083,11.1558,12.4633,12.2761,12.1107,11.9628,11.8291,11.7072,11.5952,11.4917,11.3955,11.3057,11.2214,11.1421]’;
yWT = [0,1,0,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1]’;
n = ones(size(yWT));
[b, dev, stats] = glmfit(xWT, [yWT, n], ‘binomial’, ‘Link’, ‘probit’);
wt = 5 : 0.05 : 17;
[y_fit, dylo, dyhi] = glmval(b, wt, ‘probit’, stats, ‘confidence’, 0.95);
mu_1 = -b(1) /b(2);
gradient = central_diff(y_fit, wt);
go = yWT == 1;
nogo = yWT == 0;
fig1 = figure;
plot(xWT(go), yWT(go), ‘ko’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘g’, ‘MarkerSize’, 8)
hold on
plot(xWT(nogo), yWT(nogo), ‘ks’, ‘LineWidth’, 0.5, ‘MarkerFaceColor’, ‘r’, ‘MarkerSize’, 8)
p1 = plot(wt, y_fit, ‘-m’);
p2 = plot(wt, y_fit – dylo, ‘-k’);
plot(wt, y_fit + dyhi, ‘-k’)
p3 = plot(wt, gradient/max(gradient), ‘-‘, ‘Color’, 0.6*[1,1,1]);
xlim([5, 17]);
title(‘Wu and Tian Example’)
xlabel(‘Quantile (q, in units wt)’)
ylabel(‘Probability (p)’)
grid on; grid minor
mfw_pos = get(gcf, ‘Position’); mfw_pos(3) = mfw_pos(3) * 1.6; set(gcf, ‘Position’, mfw_pos); % Make Figure Wide (60% wider)
% These are the results from the R code "Gonogo"
% plot(R_output(:,2), R_output(:,5), ‘–b’)
p4 = plot(R_output(:,2), R_output(:,6), ‘–b’);
plot(R_output(:,2), R_output(:,4), ‘–b’)
p5 = plot(R_output(:,1), R_output(:,5), ‘–r’);
plot(R_output(:,3), R_output(:,5), ‘–r’)
legend([p1 p2 p3 p4 p5],{‘Logistic Fit’,’Confidence Bounds from glmval’,’diff of Logistic Fit’,’Confidence Interval about p’,’Confidence Interval about q’}, …
‘Location’, ‘east’)
mdl = fitglm(xWT, yWT, ‘Distribution’, ‘binomial’, ‘Link’, ‘probit’);
ci = coefCI(mdl, 0.95);
figure(fig1)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,1)), ‘-c’, ‘DisplayName’, ‘from Matlab coefCI’)
plot(wt, mdl.Link.Inverse([ones(size(wt’,1),1,"like",wt’) wt’] * ci(:,2)), ‘-c’, ‘HandleVisibility’,’off’) generalized linear model, statistics, confidence limits MATLAB Answers — New Questions
Cube with different side color in RGB
I created a cube using rigid body tree (code below). now i wanted to color different side of the cube, red, gree, blue, etc… but when i use the command addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, ‘r’); for a side to be in red color, the issue is, it is giving me error of too many arguments.
cube = rigidBodyTree;
body1 = rigidBody(‘body1’);
joint1 = rigidBodyJoint(‘joint1′,’floating’);
body1.Joint = joint1;
addVisual(body1, ‘box’, [1 1 1]);
%addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, colors{i});
addBody(cube,body1,cube.BaseName);I created a cube using rigid body tree (code below). now i wanted to color different side of the cube, red, gree, blue, etc… but when i use the command addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, ‘r’); for a side to be in red color, the issue is, it is giving me error of too many arguments.
cube = rigidBodyTree;
body1 = rigidBody(‘body1’);
joint1 = rigidBodyJoint(‘joint1′,’floating’);
body1.Joint = joint1;
addVisual(body1, ‘box’, [1 1 1]);
%addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, colors{i});
addBody(cube,body1,cube.BaseName); I created a cube using rigid body tree (code below). now i wanted to color different side of the cube, red, gree, blue, etc… but when i use the command addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, ‘r’); for a side to be in red color, the issue is, it is giving me error of too many arguments.
cube = rigidBodyTree;
body1 = rigidBody(‘body1’);
joint1 = rigidBodyJoint(‘joint1′,’floating’);
body1.Joint = joint1;
addVisual(body1, ‘box’, [1 1 1]);
%addVisual(body, ‘box’, [1 1 1], ‘FaceColor’, colors{i});
addBody(cube,body1,cube.BaseName); rigid body, matlab, matlab code, addvisual, rigidbody MATLAB Answers — New Questions
[With code and PDF] How can I solve systems of PDEs with three independent variables by using the ode15s solver of the MATLAB?
[***Please find the matlab commands and PDF files***]
Now, I am trying to solve systems of PDEs (attached PDF file) by MATLAB. As you can see, this PDEs have three independent variable (z,r,t).
According to a kind suggestion by "Torsten" (http://www.mathworks.com/matlabcentral/answers/170648), I have discretized my equations about z- and r-coordinate. Then, I have developed matlab codes (PDEbook_s.m with pde_s.m, dss004,m and dss044.m) by using ode15s as follows. dss004 and dss044 (by http://www.pdecomp.net/) calculate first and second derivative values, respectively.
In this PDEs, the U and C values should decrease to zero value with increasing time. However, this values closed to a constant value. (Please run the commands in your command window) . Therefore, I think that this command cannot return correct values of the PDE systems now.
Can you find any problem in these codes?
*************PDEbook_s.m ****************************************************************************************
format short
%
% Clear previous files
clear all
clc
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
%Model parameters
zl=1;
r0=1;
%
% Grid in axial direcion
nz=10;
dz=zl/(nz-1);
for i=1:nz;
z(i)=(i-1)*dz;
end
dzs=dz^2;
%
%Grid in radial direction
nr=10;
dr=r0/(nr-1);
for j=1:nr
r(j)=(j-1)*dr;
end
drs=dr^2;
%
%Independent variables for ODE integration
tf=5;
delt=1.0E-2;
tout=[0.0:delt:(tf-delt)];
nout=(tf)/delt;
call=0;
%
%intial condition
for i=1:nz
for j=1:nr
C(i,j)=0;
U(i,j)=1;
y0((i-1)*nr+j)=C(i,j);
y0((i-1)*nr+j+nz*nr)=U(i,j);
end
end
%
%ODE integration
reltol=1.0E-4; abstol=1.0E-5;
options=odeset(‘RelTol’,reltol, ‘AbsTol’,abstol);
[t,y]=ode15s(@pde_s,tout,y0,options);
%
%1D to 2D matrices
for it=1:nout
for i=1:nz
for j=1:nr
C(it,i,j)=y(it,(i-1)*nr+j);
U(it,i,j)=y(it,(i-1)*nr+j+nz*nr);
end
end
end
%
% t-C,U plots (r=r0)
figure(1);
subplot(2,2,1)
plot(tout,C(:,1,nr),’–g’,tout,U(:,1,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=0,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,nr),’–g’,tout,U(:,5,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=4*dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,nr),’–g’,tout,U(:,8,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,nr),’–g’,tout,U(:,nz,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
%
% t-C,U plots (r=1/2r0)
figure(2);
subplot(2,2,1)
plot(tout,C(:,1,5),’–g’,tout,U(:,1,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=0,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,5),’–g’,tout,U(:,5,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=4*dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,5),’–g’,tout,U(:,8,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,5),’–g’,tout,U(:,nz,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
*****End of PDEbook_s****************************************************************************************
***********pde_s.m ****************************************************************************************
function yt=pde_s(t,y)
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
C(i,j)=y(ij);
U(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
C1d=C(i,:); %1d means 1D value
U1d=U(i,:);
%
% Ur by dss004:Ur means first derivative of U for r
Ur1d=dss004(0.0,r0,nr,U1d); %Ur means first derivative of U for r
Ur(i,:)=Ur1d;
Ur(i,1)= 0.0; %B.C.1 of Eq.(1)
Ur(i,nr)=C(i,j)-U(i,nr);%B.C.2 of Eq.(1)
%
%Urr by dss 044 with the Neumann B.D. (nl=nu=2):Urr means second derivative of U for r
Ur1d( 1)=0.0;%B.C.1 of Eq.(1)
Ur1d(nr)=C1d(j)-U1d(nr);%B.C.2 of Eq.(1)
nl=2;
nu=2;
Urr1d=dss044(0.0,r0,nr,U1d,Ur1d,nl,nu);
Urr(i,:)=Urr1d;
%
for j=1:nr
%
% (1/r)*Ur
if(j~=1)
Ur(i,j)=(1.0/r(j))*Ur(i,j);
end
if(j==1)Urr(i,j)=2.0*Urr(i,j);end
%
% Cz by dss004:Cz means first derivative of C for z
Cz1d=dss004(0.0,zl,nz,C1d);
Cz(:,j)=Cz1d;
Cz(1,j)= C1d(1);%B.C.1 of Eq.(2)
Cz(nz,j)=0.0;%B.C.2 of Eq.(2)
%
% Czz by dss044 with the Neumann B.D. (nl=nu=2):Czz means second derivative of C for z
% Czz
Cz1d(1)=C(1,nr);%B.C.1 of Eq.(2)
Cz1d(nz)=0.0;%B.C.2 of Eq.(2)
nl=2;
nu=2;
Czz1d=dss044(0.0,zl,nz,C1d,Cz1d,nl,nu);
Czz(:,j)=Czz1d;
%
% PDEs
Ut(i,j)=(1/(1+U(i,j)))*Ur(i,j)+(1/(1+U(i,j)))*Urr(i,j); %Eq.(1)
Ct(i,j)=Czz(i,j)-Cz(i,j)+(U(i,nr)-C(i,j)); % Eq.(2)
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=Ct(i,j);
yt(ij+nr*nz)=Ut(i,j);
end
end
%
% Transpose and count
yt=yt’;
ncall=ncall+1;
end
***********End of pde_s.m ****************************************************************************************
***********dss004.m ****************************************************************************************
% File: dss004.m
%
function [ux]=dss004(xl,xu,n,u)
%
% Function dss004 computes the first derivative, u , of a
% x
% variable u over the spatial domain xl le x le xu from classical
% five-point, fourth-order finite difference approximations
%
% Argument list
%
% xl Lower boundary value of x (input)
%
% xu Upper boundary value of x (input)
%
% n Number of grid points in the x domain including the
% boundary points (input)
%
% u One-dimensional array containing the values of u at
% the n grid point points for which the derivative is
% to be computed (input)
%
% ux One-dimensional array containing the numerical
% values of the derivatives of u at the n grid points
% (output)
%
% The mathematical details of the following Taylor series (or
% polynomials) are given in routine dss002.
%
% Five-point formulas
%
% (1) Left end, point i = 1
%
% 2 3 4
% a(u2 = u1 + u1 ( dx) + u1 ( dx) + u1 ( dx) + u1 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 ( dx) + u1 ( dx) + u1 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u1 + u1 (2dx) + u1 (2dx) + u1 (2dx) + u1 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (2dx) + u1 (2dx) + u1 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u1 + u1 (3dx) + u1 (3dx) + u1 (3dx) + u1 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (3dx) + u1 (3dx) + u1 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u1 + u1 (4dx) + u1 (4dx) + u1 (4dx) + u1 (4dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (4dx) + u1 (4dx) + u1 (4dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% Constants a, b, c and d are selected so that the coefficients
% of the u1 terms sum to one and the coefficients of the u1 ,
% x 2x
% u1 and u1 terms sum to zero
% 3x 4x
%
% a + 2b + 3c + 4d = 1
%
% a + 4b + 9c + 16d = 0
%
% a + 8b + 27c + 64d = 0
%
% a + 16b + 81c + 256d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u1 = (1/12dx)(-25u1 + 48u2 – 36u3 + 16u4 – 3u5) + O(dx ) (1)
% x
%
% (2) Interior point, i = 2
%
% 2 3 4
% a(u1 = u2 + u2 (-dx) + u2 (-dx) + u2 (-dx) + u2 (-dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (-dx) + u2 (-dx) + u2 (-dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u2 + u2 ( dx) + u2 ( dx) + u2 ( dx) + u2 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 ( dx) + u2 ( dx) + u2 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u2 + u2 (2dx) + u2 (2dx) + u2 (2dx) + u2 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (2dx) + u2 (2dx) + u2 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u2 + u2 (3dx) + u2 (3dx) + u2 (3dx) + u2 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (3dx) + u2 (3dx) + u2 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% -a + b + 2c + 3d = 1
%
% a + b + 4c + 9d = 0
%
% -a + b + 8c + 27d = 0
%
% a + b + 16c + 81d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u2 = (1/12dx)(-3u1 – 10u2 + 18u3 – 6u4 + u5) + O(dx ) (2)
% x
%
% (3) Interior point i, i ne 2, n-1
%
% 2 3
% a(ui-2 = ui + ui (-2dx) + ui (-2dx) + ui (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui (-2dx) + ui (-2dx) + ui (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(ui-1 = ui + ui ( -dx) + ui ( -dx) + ui ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( -dx) + ui ( -dx) + ui ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(ui+1 = ui + ui ( dx) + ui ( dx) + ui ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( dx) + ui ( dx) + ui ( dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(ui+2 = ui + ui ( 2dx) + ui ( 2dx) + ui ( 2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -2a – b + c + 2d = 1
%
% 4a + b + c + 4d = 0
%
% -8a – b + c + 8d = 0
%
% 16a + b + c + 16d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% ui = (1/12dx)(ui-2 – 8ui-1 + 0ui + 8ui+1 – ui+2) + O(dx ) (3)
% x
%
% (4) Interior point, i = n-1
%
% 2 3
% a(un-4 = un-1 + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un-1 + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un-1 + un-1 ( -dx) + un-1 (- -x) + un-1 ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( -dx) + un-1 ( -dx) + un-1 ( -dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un = un-1 + un-1 ( dx) + un-1 ( dx) + un-1 ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) + …
% 4x 4f 5x 5f 6x 6f
%
% -3a – 2b – c + d = 1
%
% 9a + 4b + c + d = 0
%
% -27a – 8b – c + d = 0
%
% 81a + 16b + c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un-1 = (1/12dx)(-un-4 + 6un-3 – 18un-2 + 10un-1 + 3un) + O(dx )
% x
% (4)
%
% (5) Right end, point i = n
%
% 2 3
% a(un-4 = un + un (-4dx) + un (-4dx) + un (-4dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-4dx) + un (-4dx) + un (-4dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un + un (-3dx) + un (-3dx) + un (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-3dx) + un (-3dx) + un (-3dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un + un (-2dx) + un (-2dx) + un (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-2dx) + un (-2dx) + un (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un-1 = un + un ( -dx) + un ( -dx) + un ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un ( -dx) + un ( -dx) + un ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -4a – 3b – 2c – d = 1
%
% 16a + 9b + 4c + d = 0
%
% -64a – 27b – 8c – d = 0
%
% 256a + 81b + 16c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un = (1/12dx)(3un-4 – 16un-3 + 36un-2 – 48un-1 + 25un) + O(dx )
% x
% (5)
%
% The weighting coefficients for equations (1) to (5) can be
% summarized as
%
% -25 48 -36 16 -3
%
% -3 -10 18 -6 1
%
% 1/12 1 -8 0 8 -1
%
% -1 6 -18 10 3
%
% 3 -16 36 -48 25
%
% which are the coefficients reported by Bickley for n = 4, m =
% 1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical
% Differentiation, Math. Gaz., vol. 25, 1941. Note – the Bickley
% coefficients have been divided by a common factor of two).
%
% Equations (1) to (5) can now be programmed to generate the
% derivative u (x) of function u(x).
% x
%
% Compute the spatial increment
dx=(xu-xl)/(n-1);
r4fdx=1./(12.*dx);
nm2=n-2;
%
% Equation (1) (note – the rhs of equations (1), (2), (3), (4)
% and (5) have been formatted so that the numerical weighting
% coefficients can be more easily associated with the Bickley
% matrix above)
ux( 1)=r4fdx*…
( -25.*u( 1) +48.*u( 2) -36.*u( 3) +16.*u( 4) -3.*u( 5));
%
% Equation (2)
ux( 2)=r4fdx*…
( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5));
%
% Equation (3)
for i=3:nm2
ux( i)=r4fdx*…
( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2));
end
%
% Equation (4)
ux(n-1)=r4fdx*…
( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n));
%
% Equation (5)
ux( n)=r4fdx*…
( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n));
***********End of dss004.m **************************************************************************************
***********dss044.m ****************************************************************************************
% File: ds044.m
%
function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%
% Function dss044 computes a fourth-order approximation of a
% second-order derivative, with or without the normal derivative
% at the boundary.
%
% Argument list
%
% xl Left value of the spatial independent variable (input)
%
% xu Right value of the spatial independent variable (input)
%
% n Number of spatial grid points, including the end
% points (input)
%
% u One-dimensional array of the dependent variable to be
% differentiated (input)
%
% ux One-dimensional array of the first derivative of u.
% The end values of ux, ux(1) and ux(n), are used in
% Neumann boundary conditions at x = xl and x = xu,
% depending on the arguments nl and nu (see the de-
% scription of nl and nu below)
%
% uxx One-dimensional array of the second derivative of u
% (output)
%
% nl Integer index for the type of boundary condition at
% x = xl (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xl
% (ux(1) is not used)
%
% 2 – Neumann boundary condition at x = xl
% (ux(1) is used)
%
% nu Integer index for the type of boundary condition at
% x = xu (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xu
% (ux(n) is not used)
%
% 2 – Neumann boundary condition at x = xu
% (ux(n) is used)
%
% The following derivation was completed by W. E. Schiesser, Depts
% of CHE and Math, Lehigh University, Bethlehem, PA 18015, USA, on
% December 15, 1986. Additional details are given in function
% dss042.
%
% ******************************************************************
%
% (1) uxx at the interior points 3, 4,…, n-2
%
% To develop a set of fourth-order correct differentiation formulas
% for the second derivative uxx, we consider first the interior
% grid points at which a symmetric formula can be used.
%
% If we consider a formula of the form
%
% a*u(i-2) + b*u(i-1) + e*u(i) + c*u(i+1) + d*u(i+2)
%
% Taylor series expansions of u(i-2), u(i-1), u(i+1) and u(i+2)
% can be substituted into this formula. We then consider the
% linear albegraic equations relating a, b, c and d which will
% retain certain terms, i.e., uxx, and drop others, e.g., uxxx,
% uxxxx and uxxxxx.
%
% Thus, for grid points 3, 4,…, n-2
%
% To retain uxx
%
% 4*a + b + c + 4*d = 2 (1)
%
% To drop uxxx
%
% -8*a – b + c + 8*d = 0 (2)
%
% To drop uxxxx
%
% 16*a + b + c + 16*d = 0 (3)
%
% To drop uxxxxx
%
% -32*a – b + c + 32*d = 0 (4)
%
% Equations (1) to (4) can be solved for a, b, c and d. If equa-
% tion (1) is added to equation (2)
%
% -4*a + 2*c + 12*d = 2 (5)
%
% If equation (1) is subtracted from equation (3)
%
% 12*a + 12*d = -2 (6)
%
% If equation (1) is added to equation (4)
%
% -28*a + 2*c + 36*d = 2 (7)
%
% Equations (5) to (7) can be solved for a, c and d. If equation
% (5) is subtracted from equation (7), and the result combined
% with equation (6)
%
% 12*a + 12*d = -2 (6)
%
% -24*a + 24*d = 0 (8)
%
% Equations (6) and (8) can be solved for a and d. From (8), a
% = d. From equation (6), a = -1/12 and d = -1/12. Then, from
% equation (5), c = 4/3, and from equation (1), b = 4/3.
%
% The final differentiation formula is then obtained as
%
% (-1/12)*u(i-2) + (4/3)*u(i-1) +
%
% (4/3)*u(i+1) + (-1/12)*u(i+2)
%
% (-1/12 + 4/3 – 1/12 + 4/3)*u(i) + uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/(12*dx**2))*(-1*u(i-2) + 16*u(i-1)
%
% – 30*u(i) + 16*u(i+1) – 1*u(i+2) (9)
%
% + O(dx**4)
%
% Note that the ux term drops out, i.e., the basic equation is
%
% -2*a – b + c + 2*d =
%
% -2*(-1/12) – (4/3) + (4/3) + 2*(-1/12) = 0
%
% Equation (9) was obtained by dropping all terms in the underlying
% Taylor series up to and including the fifth derivative, uxxxxx.
% Thus, equation (9) is exact for polynomials up to and including
% fifth order. This can be checked by substituting the functions
% 1, x, x**2, x**3, x**4 and x**5 in equation (9) and computing the
% corresponding derivatives for comparison with the known second
% derivatives. This is done for 1 merely by summing the weighting
% coefficients in equation (9), which should sum to zero, i.e.,
% -1 + 16 – 30 + 16 -1 = 0.
%
% For the remaining functions, the algebra is rather involved, but
% these functions can be checked numerically, i.e., numerical values
% of x**2, x**3, x**4 and x**5 can be substituted in equation (9)
% and the computed derivatives can be compared with the know numeri-
% cal second derivatives. This is not a proof of correctness of
% equation (9), but would likely detect any errors in equation (9).
%
% ******************************************************************
%
% (2) uxx at the interior points i = 2 and n-1
%
% For grid point 2, we consider a formula of the form
%
% a*u(i-1) + f*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) + e*u(i+4)
%
% Taylor series expansions of u(i-1), u(i+1), u(i+2), u(i+3) and
% u(i+4) when substituted into this formula give linear algebraic
% equations relating a, b, c, d and e.
%
% To drop ux
%
% -a + b + 2*c + 3*d + 4*e = 0 (10)
%
% To retain uxx
%
% a + b + 4*c + 9*d + 16*e = 2 (11)
%
% To drop uxxx
%
% -a + b + 8*c + 27*d + 64*e = 0 (12)
%
% To drop uxxxx
%
% a + b + 16*c + 81*d + 256*e = 0 (13)
%
% To drop uxxxxx
%
% -a + b + 32*c + 243*d + 1024*e = 0 (14)
%
% Equations (11), (12), (13) and (14) can be solved for a, b, c,
% d and e. If equation (10) is added to equation (11)
%
% 2*b + 6*c + 12*d +20*e = 2 (15)
%
% If equation (10) is subtracted from equation (12)
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (10) is added to equation (13)
%
% 2*b + 18*c + 84*d + 260*e = 0 (17)
%
% If equation (10) is subtracted from equation (14)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (15), (16), (17) and (18) can be solved for b, c, d
% and e.
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (15) is subtracted from equation (17)
%
% 12*c + 72*d + 240*e = -2 (19)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (16), (18) and (19) can be solved for c, d and e. If
% two times equation (16) is subtracted from equation (19),
%
% 24*d + 120*e = -2 (20)
%
% If five times equation (16) is subtracted from equation (18),
%
% 120*d + 720*e = 0 (21)
%
% Equations (20) and (21) can be solved for d and e. From (21),
% e = (-1/6)*d. Substitution in equation (20) gives d = -1/2.
% thus, e = 1/12. From equation (16), c = 7/6. From equation
% (15), b = -1/3. From equation (10), a = 5/6.
%
% The final differentiation formula is then obtained as
%
% (5/6)*u(i-1) + (-1/3)*u(i+1) + (7/6)*u(i+2) + (-1/2)*u(i+3)
%
% + (1/12)*u(i+4) = (5/6 – 1/3 + 7/6 – 1/2 + 1/12)*u(i)
%
% + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(10*u(i-1) – 15*u(i) – 4*u(i+1)
% (22)
% + 14*u(i+2) – 6*u(i+3) + 1*u(i+4)) + O(dx**4)
%
% Equation (22) will be applied at i = 2 and n-1. thus
%
% uxx(2) = (1/12*dx**2)*(10*u(1) – 15*u(2) – 4*u(3)
% (23)
% + 14*u(4) – 6*u(5) + 1*u(6)) + O(dx**4)
%
% uxx(n-1) = (1/12*dx**2)*(10*u(n) – 15*u(n-1) – 4*u(n-2)
% (24)
% + 14*u(n-3) – 6*u(n-4) + 1*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% (3) uxx at the boundary points 1 and n
%
% Finally, for grid point 1, an approximation with a Neumann bound-
% ary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*ux(i) + f*u(i)
%
% Will be used. the corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + e = 0 (25)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d = 2 (26)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d = 0 (27)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d = 0 (28)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d = 0 (29)
%
% Equations (25) to (29) can be solved for a, b, c, d and e. If
%
% Equation (26) is subtracted from equations (27), (28) and (29),
%
% 4*b + 18*c + 48*d = -2 (30)
%
% 12*b + 72*c + 240*d = -2 (31)
%
% 28*b + 234*c + 1008*d = -2 (32)
%
% Equations (30), (31) and (32) can be solved for b, c and d
%
% 18*c + 96*d = 4 (33)
%
% 108*c + 672*d = 12 (34)
%
% Equations (3) and (34) can be solved for c and d, c = 8/9, d =
% -1/8.
%
% From equation (30), b = -3. From equation (26), a = 8. From
% equation (25), e = -25/6.
%
% The final differentiation formula is then obtained as
%
% 8*u(i+1) – 3*u(i+2) + (8/9)*u(i+3) – (1/8)*u(i+4)
%
% – (25/6)*ux(i)*dx
%
% = (8 – 3 + (8/9) – (1/8))*u(i) + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*((-415/6)*u(i) + 96*u(i+1) – 36*u(i+2)
% (35)
% + (32/3)*u(i+3) – (3/2)*u(i+4) – 50*ux(i)*dx) + O(dx**4)
%
% Equation (35) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*((-415/6)*u(1) + 96*u(2) – 36*u(3)
% (36)
% + (32/3)*u(4) – (3/2)*u(5) – 50*ux(1)*dx) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*((-415/6)*u(n) + 96*u(n-1) – 36*u(n-2)
% (37)
% + (32/3)*u(n-3) – (3/2)*u(n-4) + 50*ux(n)*dx) + O(dx**4)
%
% Alternatively, for grid point 1, an approximation with a Dirichlet
% boundary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*u(i+5) + f*u(i)
%
% can be used. The corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + 5*e = 0 (38)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d + 25*e = 2 (39)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d + 125*e = 0 (40)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d + 625*e = 0 (41)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d + 3125*e = 0 (42)
%
% Equations (38), (39), (40), (41) amd (42) can be solved for a,
% b, c, d and e.
%
% 2*b + 6*c + 12*d + 20*e = 2 (43)
%
% 6*b + 24*c + 60*d + 120*e = 0 (44)
%
% 14*b + 78*c + 252*d + 620*e = 0 (45)
%
% 30*b + 240*c + 1020*d + 3120*e = 0 (46)
%
% Equations (43), (44), (45) and (46) can be solved for b, c, d
% and e
%
% 6*c + 24*d + 60*e = -6 (47)
%
% 36*c + 168*d + 480*e = -14 (48)
%
% 150*c + 840*d + 2820*e = -30 (49)
%
% Equations (47), (48) and (49) can be solved for c, d and e
%
% 24*d + 120*e = 22 (50)
%
% 240*d + 1320*e = 120 (51)
%
% From equations (50) and (51), d = 61/12, e = -5/6. From equation
% (47), c = -13. From equation (43), b = 107/6. From equation (38),
% a = -77/6.
%
% The final differentiation formula is then obtained as
%
% (-77/6)*u(i+1) + (107/6)*u(i+2) – 13*u(i+3) + (61/12)*u(i+4)
%
% – (5/6)*u(i+5) = (-77/6 + 107/6 – 13 + 61/12 – 5/6)*u(i) +
%
% uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(45*u(i) – 154*u(i+1) + 214*u(i+2)
% (52)
% – 156*u(i+3) + 61*u(i+4) – 10*u(i+5)) + O(dx**4)
%
% Equation (52) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*(45*u(1) – 154*u(2) + 214*u(3)
% (53)
% – 156*u(4) + 61*u(5) – 10*u(6)) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*(45*u(n) – 154*u(n-1) + 214*u(n-2)
% (54)
% -156*u(n-3) + 61*u(n-4) – 10*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% 1/(12*dx**2) for subsequent use
r12dxs=1./(12.0*dx^2);
%
% uxx at the left boundary
%
% Without ux (equation (53))
if nl==1
uxx(1)=r12dxs*…
( 45.0*u(1)…
-154.0*u(2)…
+214.0*u(3)…
-156.0*u(4)…
+61.0*u(5)…
-10.0*u(6));
%
% With ux (equation (36))
elseif nl==2
uxx(1)=r12dxs*…
(-415.0/6.0*u(1)…
+96.0*u(2)…
-36.0*u(3)…
+32.0/3.0*u(4)…
-3.0/2.0*u(5)…
-50.0*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux (equation (54))
if nu==1
uxx(n)=r12dxs*…
( 45.0*u(n )…
-154.0*u(n-1)…
+214.0*u(n-2)…
-156.0*u(n-3)…
+61.0*u(n-4)…
-10.0*u(n-5));
%
% With ux (equation (37))
elseif nu==2
uxx(n)=r12dxs*…
(-415.0/6.0*u(n )…
+96.0*u(n-1)…
-36.0*u(n-2)…
+32.0/3.0*u(n-3)…
-3.0/2.0*u(n-4)…
+50.0*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2 (equation (23))
uxx(2)=r12dxs*…
( 10.0*u(1)…
-15.0*u(2)…
-4.0*u(3)…
+14.0*u(4)…
-6.0*u(5)…
+1.0*u(6));
%
% i = n-1 (equation (24))
uxx(n-1)=r12dxs*…
( 10.0*u(n )…
-15.0*u(n-1)…
-4.0*u(n-2)…
+14.0*u(n-3)…
-6.0*u(n-4)…
+1.0*u(n-5));
%
% i = 3, 4,…, n-2 (equation (9))
for i=3:n-2
uxx(i)=r12dxs*…
( -1.0*u(i-2)…
+16.0*u(i-1)…
-30.0*u(i )…
+16.0*u(i+1)…
-1.0*u(i+2));
end
***********End of dss044.m **************************************************************************************[***Please find the matlab commands and PDF files***]
Now, I am trying to solve systems of PDEs (attached PDF file) by MATLAB. As you can see, this PDEs have three independent variable (z,r,t).
According to a kind suggestion by "Torsten" (http://www.mathworks.com/matlabcentral/answers/170648), I have discretized my equations about z- and r-coordinate. Then, I have developed matlab codes (PDEbook_s.m with pde_s.m, dss004,m and dss044.m) by using ode15s as follows. dss004 and dss044 (by http://www.pdecomp.net/) calculate first and second derivative values, respectively.
In this PDEs, the U and C values should decrease to zero value with increasing time. However, this values closed to a constant value. (Please run the commands in your command window) . Therefore, I think that this command cannot return correct values of the PDE systems now.
Can you find any problem in these codes?
*************PDEbook_s.m ****************************************************************************************
format short
%
% Clear previous files
clear all
clc
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
%Model parameters
zl=1;
r0=1;
%
% Grid in axial direcion
nz=10;
dz=zl/(nz-1);
for i=1:nz;
z(i)=(i-1)*dz;
end
dzs=dz^2;
%
%Grid in radial direction
nr=10;
dr=r0/(nr-1);
for j=1:nr
r(j)=(j-1)*dr;
end
drs=dr^2;
%
%Independent variables for ODE integration
tf=5;
delt=1.0E-2;
tout=[0.0:delt:(tf-delt)];
nout=(tf)/delt;
call=0;
%
%intial condition
for i=1:nz
for j=1:nr
C(i,j)=0;
U(i,j)=1;
y0((i-1)*nr+j)=C(i,j);
y0((i-1)*nr+j+nz*nr)=U(i,j);
end
end
%
%ODE integration
reltol=1.0E-4; abstol=1.0E-5;
options=odeset(‘RelTol’,reltol, ‘AbsTol’,abstol);
[t,y]=ode15s(@pde_s,tout,y0,options);
%
%1D to 2D matrices
for it=1:nout
for i=1:nz
for j=1:nr
C(it,i,j)=y(it,(i-1)*nr+j);
U(it,i,j)=y(it,(i-1)*nr+j+nz*nr);
end
end
end
%
% t-C,U plots (r=r0)
figure(1);
subplot(2,2,1)
plot(tout,C(:,1,nr),’–g’,tout,U(:,1,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=0,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,nr),’–g’,tout,U(:,5,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=4*dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,nr),’–g’,tout,U(:,8,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,nr),’–g’,tout,U(:,nz,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
%
% t-C,U plots (r=1/2r0)
figure(2);
subplot(2,2,1)
plot(tout,C(:,1,5),’–g’,tout,U(:,1,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=0,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,5),’–g’,tout,U(:,5,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=4*dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,5),’–g’,tout,U(:,8,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,5),’–g’,tout,U(:,nz,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
*****End of PDEbook_s****************************************************************************************
***********pde_s.m ****************************************************************************************
function yt=pde_s(t,y)
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
C(i,j)=y(ij);
U(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
C1d=C(i,:); %1d means 1D value
U1d=U(i,:);
%
% Ur by dss004:Ur means first derivative of U for r
Ur1d=dss004(0.0,r0,nr,U1d); %Ur means first derivative of U for r
Ur(i,:)=Ur1d;
Ur(i,1)= 0.0; %B.C.1 of Eq.(1)
Ur(i,nr)=C(i,j)-U(i,nr);%B.C.2 of Eq.(1)
%
%Urr by dss 044 with the Neumann B.D. (nl=nu=2):Urr means second derivative of U for r
Ur1d( 1)=0.0;%B.C.1 of Eq.(1)
Ur1d(nr)=C1d(j)-U1d(nr);%B.C.2 of Eq.(1)
nl=2;
nu=2;
Urr1d=dss044(0.0,r0,nr,U1d,Ur1d,nl,nu);
Urr(i,:)=Urr1d;
%
for j=1:nr
%
% (1/r)*Ur
if(j~=1)
Ur(i,j)=(1.0/r(j))*Ur(i,j);
end
if(j==1)Urr(i,j)=2.0*Urr(i,j);end
%
% Cz by dss004:Cz means first derivative of C for z
Cz1d=dss004(0.0,zl,nz,C1d);
Cz(:,j)=Cz1d;
Cz(1,j)= C1d(1);%B.C.1 of Eq.(2)
Cz(nz,j)=0.0;%B.C.2 of Eq.(2)
%
% Czz by dss044 with the Neumann B.D. (nl=nu=2):Czz means second derivative of C for z
% Czz
Cz1d(1)=C(1,nr);%B.C.1 of Eq.(2)
Cz1d(nz)=0.0;%B.C.2 of Eq.(2)
nl=2;
nu=2;
Czz1d=dss044(0.0,zl,nz,C1d,Cz1d,nl,nu);
Czz(:,j)=Czz1d;
%
% PDEs
Ut(i,j)=(1/(1+U(i,j)))*Ur(i,j)+(1/(1+U(i,j)))*Urr(i,j); %Eq.(1)
Ct(i,j)=Czz(i,j)-Cz(i,j)+(U(i,nr)-C(i,j)); % Eq.(2)
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=Ct(i,j);
yt(ij+nr*nz)=Ut(i,j);
end
end
%
% Transpose and count
yt=yt’;
ncall=ncall+1;
end
***********End of pde_s.m ****************************************************************************************
***********dss004.m ****************************************************************************************
% File: dss004.m
%
function [ux]=dss004(xl,xu,n,u)
%
% Function dss004 computes the first derivative, u , of a
% x
% variable u over the spatial domain xl le x le xu from classical
% five-point, fourth-order finite difference approximations
%
% Argument list
%
% xl Lower boundary value of x (input)
%
% xu Upper boundary value of x (input)
%
% n Number of grid points in the x domain including the
% boundary points (input)
%
% u One-dimensional array containing the values of u at
% the n grid point points for which the derivative is
% to be computed (input)
%
% ux One-dimensional array containing the numerical
% values of the derivatives of u at the n grid points
% (output)
%
% The mathematical details of the following Taylor series (or
% polynomials) are given in routine dss002.
%
% Five-point formulas
%
% (1) Left end, point i = 1
%
% 2 3 4
% a(u2 = u1 + u1 ( dx) + u1 ( dx) + u1 ( dx) + u1 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 ( dx) + u1 ( dx) + u1 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u1 + u1 (2dx) + u1 (2dx) + u1 (2dx) + u1 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (2dx) + u1 (2dx) + u1 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u1 + u1 (3dx) + u1 (3dx) + u1 (3dx) + u1 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (3dx) + u1 (3dx) + u1 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u1 + u1 (4dx) + u1 (4dx) + u1 (4dx) + u1 (4dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (4dx) + u1 (4dx) + u1 (4dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% Constants a, b, c and d are selected so that the coefficients
% of the u1 terms sum to one and the coefficients of the u1 ,
% x 2x
% u1 and u1 terms sum to zero
% 3x 4x
%
% a + 2b + 3c + 4d = 1
%
% a + 4b + 9c + 16d = 0
%
% a + 8b + 27c + 64d = 0
%
% a + 16b + 81c + 256d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u1 = (1/12dx)(-25u1 + 48u2 – 36u3 + 16u4 – 3u5) + O(dx ) (1)
% x
%
% (2) Interior point, i = 2
%
% 2 3 4
% a(u1 = u2 + u2 (-dx) + u2 (-dx) + u2 (-dx) + u2 (-dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (-dx) + u2 (-dx) + u2 (-dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u2 + u2 ( dx) + u2 ( dx) + u2 ( dx) + u2 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 ( dx) + u2 ( dx) + u2 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u2 + u2 (2dx) + u2 (2dx) + u2 (2dx) + u2 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (2dx) + u2 (2dx) + u2 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u2 + u2 (3dx) + u2 (3dx) + u2 (3dx) + u2 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (3dx) + u2 (3dx) + u2 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% -a + b + 2c + 3d = 1
%
% a + b + 4c + 9d = 0
%
% -a + b + 8c + 27d = 0
%
% a + b + 16c + 81d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u2 = (1/12dx)(-3u1 – 10u2 + 18u3 – 6u4 + u5) + O(dx ) (2)
% x
%
% (3) Interior point i, i ne 2, n-1
%
% 2 3
% a(ui-2 = ui + ui (-2dx) + ui (-2dx) + ui (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui (-2dx) + ui (-2dx) + ui (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(ui-1 = ui + ui ( -dx) + ui ( -dx) + ui ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( -dx) + ui ( -dx) + ui ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(ui+1 = ui + ui ( dx) + ui ( dx) + ui ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( dx) + ui ( dx) + ui ( dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(ui+2 = ui + ui ( 2dx) + ui ( 2dx) + ui ( 2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -2a – b + c + 2d = 1
%
% 4a + b + c + 4d = 0
%
% -8a – b + c + 8d = 0
%
% 16a + b + c + 16d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% ui = (1/12dx)(ui-2 – 8ui-1 + 0ui + 8ui+1 – ui+2) + O(dx ) (3)
% x
%
% (4) Interior point, i = n-1
%
% 2 3
% a(un-4 = un-1 + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un-1 + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un-1 + un-1 ( -dx) + un-1 (- -x) + un-1 ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( -dx) + un-1 ( -dx) + un-1 ( -dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un = un-1 + un-1 ( dx) + un-1 ( dx) + un-1 ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) + …
% 4x 4f 5x 5f 6x 6f
%
% -3a – 2b – c + d = 1
%
% 9a + 4b + c + d = 0
%
% -27a – 8b – c + d = 0
%
% 81a + 16b + c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un-1 = (1/12dx)(-un-4 + 6un-3 – 18un-2 + 10un-1 + 3un) + O(dx )
% x
% (4)
%
% (5) Right end, point i = n
%
% 2 3
% a(un-4 = un + un (-4dx) + un (-4dx) + un (-4dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-4dx) + un (-4dx) + un (-4dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un + un (-3dx) + un (-3dx) + un (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-3dx) + un (-3dx) + un (-3dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un + un (-2dx) + un (-2dx) + un (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-2dx) + un (-2dx) + un (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un-1 = un + un ( -dx) + un ( -dx) + un ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un ( -dx) + un ( -dx) + un ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -4a – 3b – 2c – d = 1
%
% 16a + 9b + 4c + d = 0
%
% -64a – 27b – 8c – d = 0
%
% 256a + 81b + 16c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un = (1/12dx)(3un-4 – 16un-3 + 36un-2 – 48un-1 + 25un) + O(dx )
% x
% (5)
%
% The weighting coefficients for equations (1) to (5) can be
% summarized as
%
% -25 48 -36 16 -3
%
% -3 -10 18 -6 1
%
% 1/12 1 -8 0 8 -1
%
% -1 6 -18 10 3
%
% 3 -16 36 -48 25
%
% which are the coefficients reported by Bickley for n = 4, m =
% 1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical
% Differentiation, Math. Gaz., vol. 25, 1941. Note – the Bickley
% coefficients have been divided by a common factor of two).
%
% Equations (1) to (5) can now be programmed to generate the
% derivative u (x) of function u(x).
% x
%
% Compute the spatial increment
dx=(xu-xl)/(n-1);
r4fdx=1./(12.*dx);
nm2=n-2;
%
% Equation (1) (note – the rhs of equations (1), (2), (3), (4)
% and (5) have been formatted so that the numerical weighting
% coefficients can be more easily associated with the Bickley
% matrix above)
ux( 1)=r4fdx*…
( -25.*u( 1) +48.*u( 2) -36.*u( 3) +16.*u( 4) -3.*u( 5));
%
% Equation (2)
ux( 2)=r4fdx*…
( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5));
%
% Equation (3)
for i=3:nm2
ux( i)=r4fdx*…
( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2));
end
%
% Equation (4)
ux(n-1)=r4fdx*…
( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n));
%
% Equation (5)
ux( n)=r4fdx*…
( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n));
***********End of dss004.m **************************************************************************************
***********dss044.m ****************************************************************************************
% File: ds044.m
%
function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%
% Function dss044 computes a fourth-order approximation of a
% second-order derivative, with or without the normal derivative
% at the boundary.
%
% Argument list
%
% xl Left value of the spatial independent variable (input)
%
% xu Right value of the spatial independent variable (input)
%
% n Number of spatial grid points, including the end
% points (input)
%
% u One-dimensional array of the dependent variable to be
% differentiated (input)
%
% ux One-dimensional array of the first derivative of u.
% The end values of ux, ux(1) and ux(n), are used in
% Neumann boundary conditions at x = xl and x = xu,
% depending on the arguments nl and nu (see the de-
% scription of nl and nu below)
%
% uxx One-dimensional array of the second derivative of u
% (output)
%
% nl Integer index for the type of boundary condition at
% x = xl (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xl
% (ux(1) is not used)
%
% 2 – Neumann boundary condition at x = xl
% (ux(1) is used)
%
% nu Integer index for the type of boundary condition at
% x = xu (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xu
% (ux(n) is not used)
%
% 2 – Neumann boundary condition at x = xu
% (ux(n) is used)
%
% The following derivation was completed by W. E. Schiesser, Depts
% of CHE and Math, Lehigh University, Bethlehem, PA 18015, USA, on
% December 15, 1986. Additional details are given in function
% dss042.
%
% ******************************************************************
%
% (1) uxx at the interior points 3, 4,…, n-2
%
% To develop a set of fourth-order correct differentiation formulas
% for the second derivative uxx, we consider first the interior
% grid points at which a symmetric formula can be used.
%
% If we consider a formula of the form
%
% a*u(i-2) + b*u(i-1) + e*u(i) + c*u(i+1) + d*u(i+2)
%
% Taylor series expansions of u(i-2), u(i-1), u(i+1) and u(i+2)
% can be substituted into this formula. We then consider the
% linear albegraic equations relating a, b, c and d which will
% retain certain terms, i.e., uxx, and drop others, e.g., uxxx,
% uxxxx and uxxxxx.
%
% Thus, for grid points 3, 4,…, n-2
%
% To retain uxx
%
% 4*a + b + c + 4*d = 2 (1)
%
% To drop uxxx
%
% -8*a – b + c + 8*d = 0 (2)
%
% To drop uxxxx
%
% 16*a + b + c + 16*d = 0 (3)
%
% To drop uxxxxx
%
% -32*a – b + c + 32*d = 0 (4)
%
% Equations (1) to (4) can be solved for a, b, c and d. If equa-
% tion (1) is added to equation (2)
%
% -4*a + 2*c + 12*d = 2 (5)
%
% If equation (1) is subtracted from equation (3)
%
% 12*a + 12*d = -2 (6)
%
% If equation (1) is added to equation (4)
%
% -28*a + 2*c + 36*d = 2 (7)
%
% Equations (5) to (7) can be solved for a, c and d. If equation
% (5) is subtracted from equation (7), and the result combined
% with equation (6)
%
% 12*a + 12*d = -2 (6)
%
% -24*a + 24*d = 0 (8)
%
% Equations (6) and (8) can be solved for a and d. From (8), a
% = d. From equation (6), a = -1/12 and d = -1/12. Then, from
% equation (5), c = 4/3, and from equation (1), b = 4/3.
%
% The final differentiation formula is then obtained as
%
% (-1/12)*u(i-2) + (4/3)*u(i-1) +
%
% (4/3)*u(i+1) + (-1/12)*u(i+2)
%
% (-1/12 + 4/3 – 1/12 + 4/3)*u(i) + uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/(12*dx**2))*(-1*u(i-2) + 16*u(i-1)
%
% – 30*u(i) + 16*u(i+1) – 1*u(i+2) (9)
%
% + O(dx**4)
%
% Note that the ux term drops out, i.e., the basic equation is
%
% -2*a – b + c + 2*d =
%
% -2*(-1/12) – (4/3) + (4/3) + 2*(-1/12) = 0
%
% Equation (9) was obtained by dropping all terms in the underlying
% Taylor series up to and including the fifth derivative, uxxxxx.
% Thus, equation (9) is exact for polynomials up to and including
% fifth order. This can be checked by substituting the functions
% 1, x, x**2, x**3, x**4 and x**5 in equation (9) and computing the
% corresponding derivatives for comparison with the known second
% derivatives. This is done for 1 merely by summing the weighting
% coefficients in equation (9), which should sum to zero, i.e.,
% -1 + 16 – 30 + 16 -1 = 0.
%
% For the remaining functions, the algebra is rather involved, but
% these functions can be checked numerically, i.e., numerical values
% of x**2, x**3, x**4 and x**5 can be substituted in equation (9)
% and the computed derivatives can be compared with the know numeri-
% cal second derivatives. This is not a proof of correctness of
% equation (9), but would likely detect any errors in equation (9).
%
% ******************************************************************
%
% (2) uxx at the interior points i = 2 and n-1
%
% For grid point 2, we consider a formula of the form
%
% a*u(i-1) + f*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) + e*u(i+4)
%
% Taylor series expansions of u(i-1), u(i+1), u(i+2), u(i+3) and
% u(i+4) when substituted into this formula give linear algebraic
% equations relating a, b, c, d and e.
%
% To drop ux
%
% -a + b + 2*c + 3*d + 4*e = 0 (10)
%
% To retain uxx
%
% a + b + 4*c + 9*d + 16*e = 2 (11)
%
% To drop uxxx
%
% -a + b + 8*c + 27*d + 64*e = 0 (12)
%
% To drop uxxxx
%
% a + b + 16*c + 81*d + 256*e = 0 (13)
%
% To drop uxxxxx
%
% -a + b + 32*c + 243*d + 1024*e = 0 (14)
%
% Equations (11), (12), (13) and (14) can be solved for a, b, c,
% d and e. If equation (10) is added to equation (11)
%
% 2*b + 6*c + 12*d +20*e = 2 (15)
%
% If equation (10) is subtracted from equation (12)
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (10) is added to equation (13)
%
% 2*b + 18*c + 84*d + 260*e = 0 (17)
%
% If equation (10) is subtracted from equation (14)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (15), (16), (17) and (18) can be solved for b, c, d
% and e.
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (15) is subtracted from equation (17)
%
% 12*c + 72*d + 240*e = -2 (19)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (16), (18) and (19) can be solved for c, d and e. If
% two times equation (16) is subtracted from equation (19),
%
% 24*d + 120*e = -2 (20)
%
% If five times equation (16) is subtracted from equation (18),
%
% 120*d + 720*e = 0 (21)
%
% Equations (20) and (21) can be solved for d and e. From (21),
% e = (-1/6)*d. Substitution in equation (20) gives d = -1/2.
% thus, e = 1/12. From equation (16), c = 7/6. From equation
% (15), b = -1/3. From equation (10), a = 5/6.
%
% The final differentiation formula is then obtained as
%
% (5/6)*u(i-1) + (-1/3)*u(i+1) + (7/6)*u(i+2) + (-1/2)*u(i+3)
%
% + (1/12)*u(i+4) = (5/6 – 1/3 + 7/6 – 1/2 + 1/12)*u(i)
%
% + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(10*u(i-1) – 15*u(i) – 4*u(i+1)
% (22)
% + 14*u(i+2) – 6*u(i+3) + 1*u(i+4)) + O(dx**4)
%
% Equation (22) will be applied at i = 2 and n-1. thus
%
% uxx(2) = (1/12*dx**2)*(10*u(1) – 15*u(2) – 4*u(3)
% (23)
% + 14*u(4) – 6*u(5) + 1*u(6)) + O(dx**4)
%
% uxx(n-1) = (1/12*dx**2)*(10*u(n) – 15*u(n-1) – 4*u(n-2)
% (24)
% + 14*u(n-3) – 6*u(n-4) + 1*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% (3) uxx at the boundary points 1 and n
%
% Finally, for grid point 1, an approximation with a Neumann bound-
% ary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*ux(i) + f*u(i)
%
% Will be used. the corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + e = 0 (25)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d = 2 (26)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d = 0 (27)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d = 0 (28)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d = 0 (29)
%
% Equations (25) to (29) can be solved for a, b, c, d and e. If
%
% Equation (26) is subtracted from equations (27), (28) and (29),
%
% 4*b + 18*c + 48*d = -2 (30)
%
% 12*b + 72*c + 240*d = -2 (31)
%
% 28*b + 234*c + 1008*d = -2 (32)
%
% Equations (30), (31) and (32) can be solved for b, c and d
%
% 18*c + 96*d = 4 (33)
%
% 108*c + 672*d = 12 (34)
%
% Equations (3) and (34) can be solved for c and d, c = 8/9, d =
% -1/8.
%
% From equation (30), b = -3. From equation (26), a = 8. From
% equation (25), e = -25/6.
%
% The final differentiation formula is then obtained as
%
% 8*u(i+1) – 3*u(i+2) + (8/9)*u(i+3) – (1/8)*u(i+4)
%
% – (25/6)*ux(i)*dx
%
% = (8 – 3 + (8/9) – (1/8))*u(i) + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*((-415/6)*u(i) + 96*u(i+1) – 36*u(i+2)
% (35)
% + (32/3)*u(i+3) – (3/2)*u(i+4) – 50*ux(i)*dx) + O(dx**4)
%
% Equation (35) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*((-415/6)*u(1) + 96*u(2) – 36*u(3)
% (36)
% + (32/3)*u(4) – (3/2)*u(5) – 50*ux(1)*dx) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*((-415/6)*u(n) + 96*u(n-1) – 36*u(n-2)
% (37)
% + (32/3)*u(n-3) – (3/2)*u(n-4) + 50*ux(n)*dx) + O(dx**4)
%
% Alternatively, for grid point 1, an approximation with a Dirichlet
% boundary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*u(i+5) + f*u(i)
%
% can be used. The corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + 5*e = 0 (38)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d + 25*e = 2 (39)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d + 125*e = 0 (40)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d + 625*e = 0 (41)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d + 3125*e = 0 (42)
%
% Equations (38), (39), (40), (41) amd (42) can be solved for a,
% b, c, d and e.
%
% 2*b + 6*c + 12*d + 20*e = 2 (43)
%
% 6*b + 24*c + 60*d + 120*e = 0 (44)
%
% 14*b + 78*c + 252*d + 620*e = 0 (45)
%
% 30*b + 240*c + 1020*d + 3120*e = 0 (46)
%
% Equations (43), (44), (45) and (46) can be solved for b, c, d
% and e
%
% 6*c + 24*d + 60*e = -6 (47)
%
% 36*c + 168*d + 480*e = -14 (48)
%
% 150*c + 840*d + 2820*e = -30 (49)
%
% Equations (47), (48) and (49) can be solved for c, d and e
%
% 24*d + 120*e = 22 (50)
%
% 240*d + 1320*e = 120 (51)
%
% From equations (50) and (51), d = 61/12, e = -5/6. From equation
% (47), c = -13. From equation (43), b = 107/6. From equation (38),
% a = -77/6.
%
% The final differentiation formula is then obtained as
%
% (-77/6)*u(i+1) + (107/6)*u(i+2) – 13*u(i+3) + (61/12)*u(i+4)
%
% – (5/6)*u(i+5) = (-77/6 + 107/6 – 13 + 61/12 – 5/6)*u(i) +
%
% uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(45*u(i) – 154*u(i+1) + 214*u(i+2)
% (52)
% – 156*u(i+3) + 61*u(i+4) – 10*u(i+5)) + O(dx**4)
%
% Equation (52) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*(45*u(1) – 154*u(2) + 214*u(3)
% (53)
% – 156*u(4) + 61*u(5) – 10*u(6)) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*(45*u(n) – 154*u(n-1) + 214*u(n-2)
% (54)
% -156*u(n-3) + 61*u(n-4) – 10*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% 1/(12*dx**2) for subsequent use
r12dxs=1./(12.0*dx^2);
%
% uxx at the left boundary
%
% Without ux (equation (53))
if nl==1
uxx(1)=r12dxs*…
( 45.0*u(1)…
-154.0*u(2)…
+214.0*u(3)…
-156.0*u(4)…
+61.0*u(5)…
-10.0*u(6));
%
% With ux (equation (36))
elseif nl==2
uxx(1)=r12dxs*…
(-415.0/6.0*u(1)…
+96.0*u(2)…
-36.0*u(3)…
+32.0/3.0*u(4)…
-3.0/2.0*u(5)…
-50.0*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux (equation (54))
if nu==1
uxx(n)=r12dxs*…
( 45.0*u(n )…
-154.0*u(n-1)…
+214.0*u(n-2)…
-156.0*u(n-3)…
+61.0*u(n-4)…
-10.0*u(n-5));
%
% With ux (equation (37))
elseif nu==2
uxx(n)=r12dxs*…
(-415.0/6.0*u(n )…
+96.0*u(n-1)…
-36.0*u(n-2)…
+32.0/3.0*u(n-3)…
-3.0/2.0*u(n-4)…
+50.0*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2 (equation (23))
uxx(2)=r12dxs*…
( 10.0*u(1)…
-15.0*u(2)…
-4.0*u(3)…
+14.0*u(4)…
-6.0*u(5)…
+1.0*u(6));
%
% i = n-1 (equation (24))
uxx(n-1)=r12dxs*…
( 10.0*u(n )…
-15.0*u(n-1)…
-4.0*u(n-2)…
+14.0*u(n-3)…
-6.0*u(n-4)…
+1.0*u(n-5));
%
% i = 3, 4,…, n-2 (equation (9))
for i=3:n-2
uxx(i)=r12dxs*…
( -1.0*u(i-2)…
+16.0*u(i-1)…
-30.0*u(i )…
+16.0*u(i+1)…
-1.0*u(i+2));
end
***********End of dss044.m ************************************************************************************** [***Please find the matlab commands and PDF files***]
Now, I am trying to solve systems of PDEs (attached PDF file) by MATLAB. As you can see, this PDEs have three independent variable (z,r,t).
According to a kind suggestion by "Torsten" (http://www.mathworks.com/matlabcentral/answers/170648), I have discretized my equations about z- and r-coordinate. Then, I have developed matlab codes (PDEbook_s.m with pde_s.m, dss004,m and dss044.m) by using ode15s as follows. dss004 and dss044 (by http://www.pdecomp.net/) calculate first and second derivative values, respectively.
In this PDEs, the U and C values should decrease to zero value with increasing time. However, this values closed to a constant value. (Please run the commands in your command window) . Therefore, I think that this command cannot return correct values of the PDE systems now.
Can you find any problem in these codes?
*************PDEbook_s.m ****************************************************************************************
format short
%
% Clear previous files
clear all
clc
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
%Model parameters
zl=1;
r0=1;
%
% Grid in axial direcion
nz=10;
dz=zl/(nz-1);
for i=1:nz;
z(i)=(i-1)*dz;
end
dzs=dz^2;
%
%Grid in radial direction
nr=10;
dr=r0/(nr-1);
for j=1:nr
r(j)=(j-1)*dr;
end
drs=dr^2;
%
%Independent variables for ODE integration
tf=5;
delt=1.0E-2;
tout=[0.0:delt:(tf-delt)];
nout=(tf)/delt;
call=0;
%
%intial condition
for i=1:nz
for j=1:nr
C(i,j)=0;
U(i,j)=1;
y0((i-1)*nr+j)=C(i,j);
y0((i-1)*nr+j+nz*nr)=U(i,j);
end
end
%
%ODE integration
reltol=1.0E-4; abstol=1.0E-5;
options=odeset(‘RelTol’,reltol, ‘AbsTol’,abstol);
[t,y]=ode15s(@pde_s,tout,y0,options);
%
%1D to 2D matrices
for it=1:nout
for i=1:nz
for j=1:nr
C(it,i,j)=y(it,(i-1)*nr+j);
U(it,i,j)=y(it,(i-1)*nr+j+nz*nr);
end
end
end
%
% t-C,U plots (r=r0)
figure(1);
subplot(2,2,1)
plot(tout,C(:,1,nr),’–g’,tout,U(:,1,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=0,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,nr),’–g’,tout,U(:,5,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=4*dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,nr),’–g’,tout,U(:,8,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,nr),’–g’,tout,U(:,nz,nr),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
%
% t-C,U plots (r=1/2r0)
figure(2);
subplot(2,2,1)
plot(tout,C(:,1,5),’–g’,tout,U(:,1,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=0,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,2)
plot(tout,C(:,5,5),’–g’,tout,U(:,5,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(=4*dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,3)
plot(tout,C(:,8,5),’–g’,tout,U(:,8,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=*7dz,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
subplot(2,2,4)
plot(tout,C(:,nz,5),’–g’,tout,U(:,nz,5),’:r’); axis([0 tf 0 1]);
title(‘Solution(z=zl,r=1/2r0)’)
xlabel(‘time’)
ylabel(‘C(green),U(red)’)
*****End of PDEbook_s****************************************************************************************
***********pde_s.m ****************************************************************************************
function yt=pde_s(t,y)
%
%Global data
global nr nz dr dz drs dzs…
r z zl r0 call tf ncall…
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
C(i,j)=y(ij);
U(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
C1d=C(i,:); %1d means 1D value
U1d=U(i,:);
%
% Ur by dss004:Ur means first derivative of U for r
Ur1d=dss004(0.0,r0,nr,U1d); %Ur means first derivative of U for r
Ur(i,:)=Ur1d;
Ur(i,1)= 0.0; %B.C.1 of Eq.(1)
Ur(i,nr)=C(i,j)-U(i,nr);%B.C.2 of Eq.(1)
%
%Urr by dss 044 with the Neumann B.D. (nl=nu=2):Urr means second derivative of U for r
Ur1d( 1)=0.0;%B.C.1 of Eq.(1)
Ur1d(nr)=C1d(j)-U1d(nr);%B.C.2 of Eq.(1)
nl=2;
nu=2;
Urr1d=dss044(0.0,r0,nr,U1d,Ur1d,nl,nu);
Urr(i,:)=Urr1d;
%
for j=1:nr
%
% (1/r)*Ur
if(j~=1)
Ur(i,j)=(1.0/r(j))*Ur(i,j);
end
if(j==1)Urr(i,j)=2.0*Urr(i,j);end
%
% Cz by dss004:Cz means first derivative of C for z
Cz1d=dss004(0.0,zl,nz,C1d);
Cz(:,j)=Cz1d;
Cz(1,j)= C1d(1);%B.C.1 of Eq.(2)
Cz(nz,j)=0.0;%B.C.2 of Eq.(2)
%
% Czz by dss044 with the Neumann B.D. (nl=nu=2):Czz means second derivative of C for z
% Czz
Cz1d(1)=C(1,nr);%B.C.1 of Eq.(2)
Cz1d(nz)=0.0;%B.C.2 of Eq.(2)
nl=2;
nu=2;
Czz1d=dss044(0.0,zl,nz,C1d,Cz1d,nl,nu);
Czz(:,j)=Czz1d;
%
% PDEs
Ut(i,j)=(1/(1+U(i,j)))*Ur(i,j)+(1/(1+U(i,j)))*Urr(i,j); %Eq.(1)
Ct(i,j)=Czz(i,j)-Cz(i,j)+(U(i,nr)-C(i,j)); % Eq.(2)
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=Ct(i,j);
yt(ij+nr*nz)=Ut(i,j);
end
end
%
% Transpose and count
yt=yt’;
ncall=ncall+1;
end
***********End of pde_s.m ****************************************************************************************
***********dss004.m ****************************************************************************************
% File: dss004.m
%
function [ux]=dss004(xl,xu,n,u)
%
% Function dss004 computes the first derivative, u , of a
% x
% variable u over the spatial domain xl le x le xu from classical
% five-point, fourth-order finite difference approximations
%
% Argument list
%
% xl Lower boundary value of x (input)
%
% xu Upper boundary value of x (input)
%
% n Number of grid points in the x domain including the
% boundary points (input)
%
% u One-dimensional array containing the values of u at
% the n grid point points for which the derivative is
% to be computed (input)
%
% ux One-dimensional array containing the numerical
% values of the derivatives of u at the n grid points
% (output)
%
% The mathematical details of the following Taylor series (or
% polynomials) are given in routine dss002.
%
% Five-point formulas
%
% (1) Left end, point i = 1
%
% 2 3 4
% a(u2 = u1 + u1 ( dx) + u1 ( dx) + u1 ( dx) + u1 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 ( dx) + u1 ( dx) + u1 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u1 + u1 (2dx) + u1 (2dx) + u1 (2dx) + u1 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (2dx) + u1 (2dx) + u1 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u1 + u1 (3dx) + u1 (3dx) + u1 (3dx) + u1 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (3dx) + u1 (3dx) + u1 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u1 + u1 (4dx) + u1 (4dx) + u1 (4dx) + u1 (4dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (4dx) + u1 (4dx) + u1 (4dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% Constants a, b, c and d are selected so that the coefficients
% of the u1 terms sum to one and the coefficients of the u1 ,
% x 2x
% u1 and u1 terms sum to zero
% 3x 4x
%
% a + 2b + 3c + 4d = 1
%
% a + 4b + 9c + 16d = 0
%
% a + 8b + 27c + 64d = 0
%
% a + 16b + 81c + 256d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u1 = (1/12dx)(-25u1 + 48u2 – 36u3 + 16u4 – 3u5) + O(dx ) (1)
% x
%
% (2) Interior point, i = 2
%
% 2 3 4
% a(u1 = u2 + u2 (-dx) + u2 (-dx) + u2 (-dx) + u2 (-dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (-dx) + u2 (-dx) + u2 (-dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u2 + u2 ( dx) + u2 ( dx) + u2 ( dx) + u2 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 ( dx) + u2 ( dx) + u2 ( dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u2 + u2 (2dx) + u2 (2dx) + u2 (2dx) + u2 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (2dx) + u2 (2dx) + u2 (2dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u2 + u2 (3dx) + u2 (3dx) + u2 (3dx) + u2 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (3dx) + u2 (3dx) + u2 (3dx) + …)
% 5x 5f 6x 6f 7x 7f
%
% -a + b + 2c + 3d = 1
%
% a + b + 4c + 9d = 0
%
% -a + b + 8c + 27d = 0
%
% a + b + 16c + 81d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u2 = (1/12dx)(-3u1 – 10u2 + 18u3 – 6u4 + u5) + O(dx ) (2)
% x
%
% (3) Interior point i, i ne 2, n-1
%
% 2 3
% a(ui-2 = ui + ui (-2dx) + ui (-2dx) + ui (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui (-2dx) + ui (-2dx) + ui (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(ui-1 = ui + ui ( -dx) + ui ( -dx) + ui ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( -dx) + ui ( -dx) + ui ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(ui+1 = ui + ui ( dx) + ui ( dx) + ui ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( dx) + ui ( dx) + ui ( dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(ui+2 = ui + ui ( 2dx) + ui ( 2dx) + ui ( 2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -2a – b + c + 2d = 1
%
% 4a + b + c + 4d = 0
%
% -8a – b + c + 8d = 0
%
% 16a + b + c + 16d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% ui = (1/12dx)(ui-2 – 8ui-1 + 0ui + 8ui+1 – ui+2) + O(dx ) (3)
% x
%
% (4) Interior point, i = n-1
%
% 2 3
% a(un-4 = un-1 + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un-1 + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un-1 + un-1 ( -dx) + un-1 (- -x) + un-1 ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( -dx) + un-1 ( -dx) + un-1 ( -dx) + …
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un = un-1 + un-1 ( dx) + un-1 ( dx) + un-1 ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) + …
% 4x 4f 5x 5f 6x 6f
%
% -3a – 2b – c + d = 1
%
% 9a + 4b + c + d = 0
%
% -27a – 8b – c + d = 0
%
% 81a + 16b + c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un-1 = (1/12dx)(-un-4 + 6un-3 – 18un-2 + 10un-1 + 3un) + O(dx )
% x
% (4)
%
% (5) Right end, point i = n
%
% 2 3
% a(un-4 = un + un (-4dx) + un (-4dx) + un (-4dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-4dx) + un (-4dx) + un (-4dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un + un (-3dx) + un (-3dx) + un (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-3dx) + un (-3dx) + un (-3dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un + un (-2dx) + un (-2dx) + un (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-2dx) + un (-2dx) + un (-2dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un-1 = un + un ( -dx) + un ( -dx) + un ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un ( -dx) + un ( -dx) + un ( -dx) + …)
% 4x 4f 5x 5f 6x 6f
%
% -4a – 3b – 2c – d = 1
%
% 16a + 9b + 4c + d = 0
%
% -64a – 27b – 8c – d = 0
%
% 256a + 81b + 16c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un = (1/12dx)(3un-4 – 16un-3 + 36un-2 – 48un-1 + 25un) + O(dx )
% x
% (5)
%
% The weighting coefficients for equations (1) to (5) can be
% summarized as
%
% -25 48 -36 16 -3
%
% -3 -10 18 -6 1
%
% 1/12 1 -8 0 8 -1
%
% -1 6 -18 10 3
%
% 3 -16 36 -48 25
%
% which are the coefficients reported by Bickley for n = 4, m =
% 1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical
% Differentiation, Math. Gaz., vol. 25, 1941. Note – the Bickley
% coefficients have been divided by a common factor of two).
%
% Equations (1) to (5) can now be programmed to generate the
% derivative u (x) of function u(x).
% x
%
% Compute the spatial increment
dx=(xu-xl)/(n-1);
r4fdx=1./(12.*dx);
nm2=n-2;
%
% Equation (1) (note – the rhs of equations (1), (2), (3), (4)
% and (5) have been formatted so that the numerical weighting
% coefficients can be more easily associated with the Bickley
% matrix above)
ux( 1)=r4fdx*…
( -25.*u( 1) +48.*u( 2) -36.*u( 3) +16.*u( 4) -3.*u( 5));
%
% Equation (2)
ux( 2)=r4fdx*…
( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5));
%
% Equation (3)
for i=3:nm2
ux( i)=r4fdx*…
( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2));
end
%
% Equation (4)
ux(n-1)=r4fdx*…
( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n));
%
% Equation (5)
ux( n)=r4fdx*…
( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n));
***********End of dss004.m **************************************************************************************
***********dss044.m ****************************************************************************************
% File: ds044.m
%
function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%
% Function dss044 computes a fourth-order approximation of a
% second-order derivative, with or without the normal derivative
% at the boundary.
%
% Argument list
%
% xl Left value of the spatial independent variable (input)
%
% xu Right value of the spatial independent variable (input)
%
% n Number of spatial grid points, including the end
% points (input)
%
% u One-dimensional array of the dependent variable to be
% differentiated (input)
%
% ux One-dimensional array of the first derivative of u.
% The end values of ux, ux(1) and ux(n), are used in
% Neumann boundary conditions at x = xl and x = xu,
% depending on the arguments nl and nu (see the de-
% scription of nl and nu below)
%
% uxx One-dimensional array of the second derivative of u
% (output)
%
% nl Integer index for the type of boundary condition at
% x = xl (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xl
% (ux(1) is not used)
%
% 2 – Neumann boundary condition at x = xl
% (ux(1) is used)
%
% nu Integer index for the type of boundary condition at
% x = xu (input). The allowable values are
%
% 1 – Dirichlet boundary condition at x = xu
% (ux(n) is not used)
%
% 2 – Neumann boundary condition at x = xu
% (ux(n) is used)
%
% The following derivation was completed by W. E. Schiesser, Depts
% of CHE and Math, Lehigh University, Bethlehem, PA 18015, USA, on
% December 15, 1986. Additional details are given in function
% dss042.
%
% ******************************************************************
%
% (1) uxx at the interior points 3, 4,…, n-2
%
% To develop a set of fourth-order correct differentiation formulas
% for the second derivative uxx, we consider first the interior
% grid points at which a symmetric formula can be used.
%
% If we consider a formula of the form
%
% a*u(i-2) + b*u(i-1) + e*u(i) + c*u(i+1) + d*u(i+2)
%
% Taylor series expansions of u(i-2), u(i-1), u(i+1) and u(i+2)
% can be substituted into this formula. We then consider the
% linear albegraic equations relating a, b, c and d which will
% retain certain terms, i.e., uxx, and drop others, e.g., uxxx,
% uxxxx and uxxxxx.
%
% Thus, for grid points 3, 4,…, n-2
%
% To retain uxx
%
% 4*a + b + c + 4*d = 2 (1)
%
% To drop uxxx
%
% -8*a – b + c + 8*d = 0 (2)
%
% To drop uxxxx
%
% 16*a + b + c + 16*d = 0 (3)
%
% To drop uxxxxx
%
% -32*a – b + c + 32*d = 0 (4)
%
% Equations (1) to (4) can be solved for a, b, c and d. If equa-
% tion (1) is added to equation (2)
%
% -4*a + 2*c + 12*d = 2 (5)
%
% If equation (1) is subtracted from equation (3)
%
% 12*a + 12*d = -2 (6)
%
% If equation (1) is added to equation (4)
%
% -28*a + 2*c + 36*d = 2 (7)
%
% Equations (5) to (7) can be solved for a, c and d. If equation
% (5) is subtracted from equation (7), and the result combined
% with equation (6)
%
% 12*a + 12*d = -2 (6)
%
% -24*a + 24*d = 0 (8)
%
% Equations (6) and (8) can be solved for a and d. From (8), a
% = d. From equation (6), a = -1/12 and d = -1/12. Then, from
% equation (5), c = 4/3, and from equation (1), b = 4/3.
%
% The final differentiation formula is then obtained as
%
% (-1/12)*u(i-2) + (4/3)*u(i-1) +
%
% (4/3)*u(i+1) + (-1/12)*u(i+2)
%
% (-1/12 + 4/3 – 1/12 + 4/3)*u(i) + uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/(12*dx**2))*(-1*u(i-2) + 16*u(i-1)
%
% – 30*u(i) + 16*u(i+1) – 1*u(i+2) (9)
%
% + O(dx**4)
%
% Note that the ux term drops out, i.e., the basic equation is
%
% -2*a – b + c + 2*d =
%
% -2*(-1/12) – (4/3) + (4/3) + 2*(-1/12) = 0
%
% Equation (9) was obtained by dropping all terms in the underlying
% Taylor series up to and including the fifth derivative, uxxxxx.
% Thus, equation (9) is exact for polynomials up to and including
% fifth order. This can be checked by substituting the functions
% 1, x, x**2, x**3, x**4 and x**5 in equation (9) and computing the
% corresponding derivatives for comparison with the known second
% derivatives. This is done for 1 merely by summing the weighting
% coefficients in equation (9), which should sum to zero, i.e.,
% -1 + 16 – 30 + 16 -1 = 0.
%
% For the remaining functions, the algebra is rather involved, but
% these functions can be checked numerically, i.e., numerical values
% of x**2, x**3, x**4 and x**5 can be substituted in equation (9)
% and the computed derivatives can be compared with the know numeri-
% cal second derivatives. This is not a proof of correctness of
% equation (9), but would likely detect any errors in equation (9).
%
% ******************************************************************
%
% (2) uxx at the interior points i = 2 and n-1
%
% For grid point 2, we consider a formula of the form
%
% a*u(i-1) + f*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) + e*u(i+4)
%
% Taylor series expansions of u(i-1), u(i+1), u(i+2), u(i+3) and
% u(i+4) when substituted into this formula give linear algebraic
% equations relating a, b, c, d and e.
%
% To drop ux
%
% -a + b + 2*c + 3*d + 4*e = 0 (10)
%
% To retain uxx
%
% a + b + 4*c + 9*d + 16*e = 2 (11)
%
% To drop uxxx
%
% -a + b + 8*c + 27*d + 64*e = 0 (12)
%
% To drop uxxxx
%
% a + b + 16*c + 81*d + 256*e = 0 (13)
%
% To drop uxxxxx
%
% -a + b + 32*c + 243*d + 1024*e = 0 (14)
%
% Equations (11), (12), (13) and (14) can be solved for a, b, c,
% d and e. If equation (10) is added to equation (11)
%
% 2*b + 6*c + 12*d +20*e = 2 (15)
%
% If equation (10) is subtracted from equation (12)
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (10) is added to equation (13)
%
% 2*b + 18*c + 84*d + 260*e = 0 (17)
%
% If equation (10) is subtracted from equation (14)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (15), (16), (17) and (18) can be solved for b, c, d
% and e.
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (15) is subtracted from equation (17)
%
% 12*c + 72*d + 240*e = -2 (19)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (16), (18) and (19) can be solved for c, d and e. If
% two times equation (16) is subtracted from equation (19),
%
% 24*d + 120*e = -2 (20)
%
% If five times equation (16) is subtracted from equation (18),
%
% 120*d + 720*e = 0 (21)
%
% Equations (20) and (21) can be solved for d and e. From (21),
% e = (-1/6)*d. Substitution in equation (20) gives d = -1/2.
% thus, e = 1/12. From equation (16), c = 7/6. From equation
% (15), b = -1/3. From equation (10), a = 5/6.
%
% The final differentiation formula is then obtained as
%
% (5/6)*u(i-1) + (-1/3)*u(i+1) + (7/6)*u(i+2) + (-1/2)*u(i+3)
%
% + (1/12)*u(i+4) = (5/6 – 1/3 + 7/6 – 1/2 + 1/12)*u(i)
%
% + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(10*u(i-1) – 15*u(i) – 4*u(i+1)
% (22)
% + 14*u(i+2) – 6*u(i+3) + 1*u(i+4)) + O(dx**4)
%
% Equation (22) will be applied at i = 2 and n-1. thus
%
% uxx(2) = (1/12*dx**2)*(10*u(1) – 15*u(2) – 4*u(3)
% (23)
% + 14*u(4) – 6*u(5) + 1*u(6)) + O(dx**4)
%
% uxx(n-1) = (1/12*dx**2)*(10*u(n) – 15*u(n-1) – 4*u(n-2)
% (24)
% + 14*u(n-3) – 6*u(n-4) + 1*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% (3) uxx at the boundary points 1 and n
%
% Finally, for grid point 1, an approximation with a Neumann bound-
% ary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*ux(i) + f*u(i)
%
% Will be used. the corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + e = 0 (25)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d = 2 (26)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d = 0 (27)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d = 0 (28)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d = 0 (29)
%
% Equations (25) to (29) can be solved for a, b, c, d and e. If
%
% Equation (26) is subtracted from equations (27), (28) and (29),
%
% 4*b + 18*c + 48*d = -2 (30)
%
% 12*b + 72*c + 240*d = -2 (31)
%
% 28*b + 234*c + 1008*d = -2 (32)
%
% Equations (30), (31) and (32) can be solved for b, c and d
%
% 18*c + 96*d = 4 (33)
%
% 108*c + 672*d = 12 (34)
%
% Equations (3) and (34) can be solved for c and d, c = 8/9, d =
% -1/8.
%
% From equation (30), b = -3. From equation (26), a = 8. From
% equation (25), e = -25/6.
%
% The final differentiation formula is then obtained as
%
% 8*u(i+1) – 3*u(i+2) + (8/9)*u(i+3) – (1/8)*u(i+4)
%
% – (25/6)*ux(i)*dx
%
% = (8 – 3 + (8/9) – (1/8))*u(i) + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*((-415/6)*u(i) + 96*u(i+1) – 36*u(i+2)
% (35)
% + (32/3)*u(i+3) – (3/2)*u(i+4) – 50*ux(i)*dx) + O(dx**4)
%
% Equation (35) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*((-415/6)*u(1) + 96*u(2) – 36*u(3)
% (36)
% + (32/3)*u(4) – (3/2)*u(5) – 50*ux(1)*dx) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*((-415/6)*u(n) + 96*u(n-1) – 36*u(n-2)
% (37)
% + (32/3)*u(n-3) – (3/2)*u(n-4) + 50*ux(n)*dx) + O(dx**4)
%
% Alternatively, for grid point 1, an approximation with a Dirichlet
% boundary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*u(i+5) + f*u(i)
%
% can be used. The corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + 5*e = 0 (38)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d + 25*e = 2 (39)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d + 125*e = 0 (40)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d + 625*e = 0 (41)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d + 3125*e = 0 (42)
%
% Equations (38), (39), (40), (41) amd (42) can be solved for a,
% b, c, d and e.
%
% 2*b + 6*c + 12*d + 20*e = 2 (43)
%
% 6*b + 24*c + 60*d + 120*e = 0 (44)
%
% 14*b + 78*c + 252*d + 620*e = 0 (45)
%
% 30*b + 240*c + 1020*d + 3120*e = 0 (46)
%
% Equations (43), (44), (45) and (46) can be solved for b, c, d
% and e
%
% 6*c + 24*d + 60*e = -6 (47)
%
% 36*c + 168*d + 480*e = -14 (48)
%
% 150*c + 840*d + 2820*e = -30 (49)
%
% Equations (47), (48) and (49) can be solved for c, d and e
%
% 24*d + 120*e = 22 (50)
%
% 240*d + 1320*e = 120 (51)
%
% From equations (50) and (51), d = 61/12, e = -5/6. From equation
% (47), c = -13. From equation (43), b = 107/6. From equation (38),
% a = -77/6.
%
% The final differentiation formula is then obtained as
%
% (-77/6)*u(i+1) + (107/6)*u(i+2) – 13*u(i+3) + (61/12)*u(i+4)
%
% – (5/6)*u(i+5) = (-77/6 + 107/6 – 13 + 61/12 – 5/6)*u(i) +
%
% uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(45*u(i) – 154*u(i+1) + 214*u(i+2)
% (52)
% – 156*u(i+3) + 61*u(i+4) – 10*u(i+5)) + O(dx**4)
%
% Equation (52) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*(45*u(1) – 154*u(2) + 214*u(3)
% (53)
% – 156*u(4) + 61*u(5) – 10*u(6)) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*(45*u(n) – 154*u(n-1) + 214*u(n-2)
% (54)
% -156*u(n-3) + 61*u(n-4) – 10*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% 1/(12*dx**2) for subsequent use
r12dxs=1./(12.0*dx^2);
%
% uxx at the left boundary
%
% Without ux (equation (53))
if nl==1
uxx(1)=r12dxs*…
( 45.0*u(1)…
-154.0*u(2)…
+214.0*u(3)…
-156.0*u(4)…
+61.0*u(5)…
-10.0*u(6));
%
% With ux (equation (36))
elseif nl==2
uxx(1)=r12dxs*…
(-415.0/6.0*u(1)…
+96.0*u(2)…
-36.0*u(3)…
+32.0/3.0*u(4)…
-3.0/2.0*u(5)…
-50.0*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux (equation (54))
if nu==1
uxx(n)=r12dxs*…
( 45.0*u(n )…
-154.0*u(n-1)…
+214.0*u(n-2)…
-156.0*u(n-3)…
+61.0*u(n-4)…
-10.0*u(n-5));
%
% With ux (equation (37))
elseif nu==2
uxx(n)=r12dxs*…
(-415.0/6.0*u(n )…
+96.0*u(n-1)…
-36.0*u(n-2)…
+32.0/3.0*u(n-3)…
-3.0/2.0*u(n-4)…
+50.0*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2 (equation (23))
uxx(2)=r12dxs*…
( 10.0*u(1)…
-15.0*u(2)…
-4.0*u(3)…
+14.0*u(4)…
-6.0*u(5)…
+1.0*u(6));
%
% i = n-1 (equation (24))
uxx(n-1)=r12dxs*…
( 10.0*u(n )…
-15.0*u(n-1)…
-4.0*u(n-2)…
+14.0*u(n-3)…
-6.0*u(n-4)…
+1.0*u(n-5));
%
% i = 3, 4,…, n-2 (equation (9))
for i=3:n-2
uxx(i)=r12dxs*…
( -1.0*u(i-2)…
+16.0*u(i-1)…
-30.0*u(i )…
+16.0*u(i+1)…
-1.0*u(i+2));
end
***********End of dss044.m ************************************************************************************** pde, ode15s, systems of pde, partial differential equation MATLAB Answers — New Questions
I keep getting a “cannot convert double value “…” to a handle. Does anybody know how I can fix my code in order for this to go away
DOES ANYBODY KNOW HOW TO FIX THIS? I keep getting a “cannot convert double value “…” to a handle. Does anybody know how I can fix my code in order for this to go away?DOES ANYBODY KNOW HOW TO FIX THIS? I keep getting a “cannot convert double value “…” to a handle. Does anybody know how I can fix my code in order for this to go away? DOES ANYBODY KNOW HOW TO FIX THIS? I keep getting a “cannot convert double value “…” to a handle. Does anybody know how I can fix my code in order for this to go away? final project, help me, ??? MATLAB Answers — New Questions
Code generation with VNT in release mode
Hi,
I am encountering an issue when I use the embedded coder with the vehicle-network-toolbox and build in release mode. When built in release mode the generated code does not seem to actually pack any CAN messages. When built in debug mode everything works as expected.
We are using the CAN PACK and CAN UNPACK functions with a .dbc file to specify the messages. We specify the message data as a bus, which we are transferring to/from our CAN driver using hand written C++ code (not using the simulink CAN Transmit / CAN Receive blocks).Hi,
I am encountering an issue when I use the embedded coder with the vehicle-network-toolbox and build in release mode. When built in release mode the generated code does not seem to actually pack any CAN messages. When built in debug mode everything works as expected.
We are using the CAN PACK and CAN UNPACK functions with a .dbc file to specify the messages. We specify the message data as a bus, which we are transferring to/from our CAN driver using hand written C++ code (not using the simulink CAN Transmit / CAN Receive blocks). Hi,
I am encountering an issue when I use the embedded coder with the vehicle-network-toolbox and build in release mode. When built in release mode the generated code does not seem to actually pack any CAN messages. When built in debug mode everything works as expected.
We are using the CAN PACK and CAN UNPACK functions with a .dbc file to specify the messages. We specify the message data as a bus, which we are transferring to/from our CAN driver using hand written C++ code (not using the simulink CAN Transmit / CAN Receive blocks). vehicle network toolbox, can, release mode, debug mode MATLAB Answers — New Questions
Error in AWGN in MATLAB
How do I correct the following error in MATLAB? MATLAB code attached.
Error using +
Integers can only be combined with integers of the same class, or scalar doubles.
Error in awgn (line 157)
y = sig + noise;
Error in Rev2 (line 20)
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray(:), 8, ‘left-msb’).’, [], 1);
% BPSK modulation
Eb_No_low = 0;
Eb_No_high = 4;
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_low);
title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(lenna_bits, SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_high);
title(‘Received Image (4 dB SNR)’);How do I correct the following error in MATLAB? MATLAB code attached.
Error using +
Integers can only be combined with integers of the same class, or scalar doubles.
Error in awgn (line 157)
y = sig + noise;
Error in Rev2 (line 20)
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray(:), 8, ‘left-msb’).’, [], 1);
% BPSK modulation
Eb_No_low = 0;
Eb_No_high = 4;
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_low);
title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(lenna_bits, SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_high);
title(‘Received Image (4 dB SNR)’); How do I correct the following error in MATLAB? MATLAB code attached.
Error using +
Integers can only be combined with integers of the same class, or scalar doubles.
Error in awgn (line 157)
y = sig + noise;
Error in Rev2 (line 20)
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
clc;
clear all;
% Load the lenna image
lenna = imread(‘lenna.png’);
% Convert image to grayscale
lenna_gray = rgb2gray(lenna);
% Convert pixel values to bits
lenna_bits = reshape(de2bi(lenna_gray(:), 8, ‘left-msb’).’, [], 1);
% BPSK modulation
Eb_No_low = 0;
Eb_No_high = 4;
SNR_low = 10^(Eb_No_low/10);
SNR_high = 10^(Eb_No_high/10);
% Transmit and receive at low SNR
received_low = awgn(lenna_bits, SNR_low, ‘measured’);
% Demodulation
decoded_low = received_low < 0;
% Reshape decoded bits to original image size
decoded_image_low = reshape(decoded_low, size(lenna_gray));
% Plot original and received image at low SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_low);
title(‘Received Image (0 dB SNR)’);
% Transmit and receive at high SNR
received_high = awgn(lenna_bits, SNR_high, ‘measured’);
% Demodulation
decoded_high = received_high < 0;
% Reshape decoded bits to original image size
decoded_image_high = reshape(decoded_high, size(lenna_gray));
% Plot original and received image at high SNR
figure;
subplot(1,2,1);
imshow(lenna_gray);
title(‘Original Image’);
subplot(1,2,2);
imshow(decoded_image_high);
title(‘Received Image (4 dB SNR)’); awgn error MATLAB Answers — New Questions
Time vs displacement plot of a Transfer function
I am trying to plot time vs displacement plot of a transfer function. It is a simple case of spherical body in a liquid hence excitation of a spherical body is being resisted by the (elastic and viscous properties of the liquid) Hence F excitation= F inertia + F (elastic) + F viscous.
The final transfer function is of the form below where Delta Q is the relative displacement between the body and the liquid where as Vm is the excitation velocity.
I have used this transfer function to find eigen frequency for this system. Now i wants to check if the system undergoes resonance at that eigen frequency. How to do it using matlabI am trying to plot time vs displacement plot of a transfer function. It is a simple case of spherical body in a liquid hence excitation of a spherical body is being resisted by the (elastic and viscous properties of the liquid) Hence F excitation= F inertia + F (elastic) + F viscous.
The final transfer function is of the form below where Delta Q is the relative displacement between the body and the liquid where as Vm is the excitation velocity.
I have used this transfer function to find eigen frequency for this system. Now i wants to check if the system undergoes resonance at that eigen frequency. How to do it using matlab I am trying to plot time vs displacement plot of a transfer function. It is a simple case of spherical body in a liquid hence excitation of a spherical body is being resisted by the (elastic and viscous properties of the liquid) Hence F excitation= F inertia + F (elastic) + F viscous.
The final transfer function is of the form below where Delta Q is the relative displacement between the body and the liquid where as Vm is the excitation velocity.
I have used this transfer function to find eigen frequency for this system. Now i wants to check if the system undergoes resonance at that eigen frequency. How to do it using matlab matlab, plot MATLAB Answers — New Questions
Can I use the Property Inspector to edit axis labels?
I am trying to edit the Y axis labels in a boxplot from a Kruskal Wallis test. Is this possible in the Property Inspector? I have tried altering the labels and tick marks in the Property Inspector window although the changes have not been applied to the figure.I am trying to edit the Y axis labels in a boxplot from a Kruskal Wallis test. Is this possible in the Property Inspector? I have tried altering the labels and tick marks in the Property Inspector window although the changes have not been applied to the figure. I am trying to edit the Y axis labels in a boxplot from a Kruskal Wallis test. Is this possible in the Property Inspector? I have tried altering the labels and tick marks in the Property Inspector window although the changes have not been applied to the figure. property inspector, axis labels, boxplot MATLAB Answers — New Questions
Transfer Learning shows no enhancements in runtime
When configuring my network for transfer learning, to "freeze" the layers I am setting the ‘WeightLearnRateFactor’ and ‘BiasLearnRateFactor’ to 0, which from my understanding is Matlab’s intended implementation of transfer learning. However, given the same number of iterations and epochs, training my network takes the same amount of time as it would during standard training. Are there any additional steps that must be taken to implement transfer learning in a custom neural network.When configuring my network for transfer learning, to "freeze" the layers I am setting the ‘WeightLearnRateFactor’ and ‘BiasLearnRateFactor’ to 0, which from my understanding is Matlab’s intended implementation of transfer learning. However, given the same number of iterations and epochs, training my network takes the same amount of time as it would during standard training. Are there any additional steps that must be taken to implement transfer learning in a custom neural network. When configuring my network for transfer learning, to "freeze" the layers I am setting the ‘WeightLearnRateFactor’ and ‘BiasLearnRateFactor’ to 0, which from my understanding is Matlab’s intended implementation of transfer learning. However, given the same number of iterations and epochs, training my network takes the same amount of time as it would during standard training. Are there any additional steps that must be taken to implement transfer learning in a custom neural network. machine learning, neural network MATLAB Answers — New Questions
How to Convert a Scalar Struct with Vector Fields to a Vector Struct with Scalar Fields?
I have a scalar struct S with each field a row vector of the same number of elements.
S.x = [1 2];
S.y = [10 20];
I want create a struct array, A, with same number of elements as in the fields in S, with the same field names as S, where each field value in each element of A is correspondingly indexed from the original field in S. I’m not sure I’m explaining this well.
So far I have this, which works:
C = [fieldnames(s) , arrayfun(@(c) mat2cell(c{1},1,ones(1,numel(c{1}))),struct2cell(s),’Uni’,false)].’;
A = struct(C{:})
[A.x]
[A.y]
Open to suggestions on better (e.g., easier to understand) alternatives.I have a scalar struct S with each field a row vector of the same number of elements.
S.x = [1 2];
S.y = [10 20];
I want create a struct array, A, with same number of elements as in the fields in S, with the same field names as S, where each field value in each element of A is correspondingly indexed from the original field in S. I’m not sure I’m explaining this well.
So far I have this, which works:
C = [fieldnames(s) , arrayfun(@(c) mat2cell(c{1},1,ones(1,numel(c{1}))),struct2cell(s),’Uni’,false)].’;
A = struct(C{:})
[A.x]
[A.y]
Open to suggestions on better (e.g., easier to understand) alternatives. I have a scalar struct S with each field a row vector of the same number of elements.
S.x = [1 2];
S.y = [10 20];
I want create a struct array, A, with same number of elements as in the fields in S, with the same field names as S, where each field value in each element of A is correspondingly indexed from the original field in S. I’m not sure I’m explaining this well.
So far I have this, which works:
C = [fieldnames(s) , arrayfun(@(c) mat2cell(c{1},1,ones(1,numel(c{1}))),struct2cell(s),’Uni’,false)].’;
A = struct(C{:})
[A.x]
[A.y]
Open to suggestions on better (e.g., easier to understand) alternatives. struct MATLAB Answers — New Questions
Why matlab doesn’t recognize the command ”vrinstall -install”?
I’m trying to install the V-Realm builder. I have already installed the 3D world editor. When I run the script: vrinstall -install, matlab says: Unrecognized function or variable ‘vrinstall’. How can i solve this problem?I’m trying to install the V-Realm builder. I have already installed the 3D world editor. When I run the script: vrinstall -install, matlab says: Unrecognized function or variable ‘vrinstall’. How can i solve this problem? I’m trying to install the V-Realm builder. I have already installed the 3D world editor. When I run the script: vrinstall -install, matlab says: Unrecognized function or variable ‘vrinstall’. How can i solve this problem? vrinstall MATLAB Answers — New Questions
Multiply tall array with array
I have a tall array, that contains audio samples created from a dataset like this:
dataS = tall(fileDatastore(fullfile(appData.OutputDirectoryPath, ‘tmp_*.mat’), …
‘ReadFcn’, @(fname) getfield(load(fname), ‘resampled’), …
‘UniformRead’, true));
This data is large, but its element count is known. I need to multiply this data with a signal:
t = ((0:N-1)’/fsn);
carrier = exp(1j*2*pi*carrierFreq*t);
ccy = dataS.*carrier;
The last line fails, with the following error:
Error using .*
Incompatible non-scalar tall array arguments. Each of the tall arrays must be the same size in the first dimension, must be derived from a single tall array, and must not have been indexed differently in the first
dimension (indexing operations include functions such as VERTCAT, SPLITAPPLY, SORT, CELL2MAT, SYNCHRONIZE, RETIME and so on).
I assume the issue is that it cannot guarantee that the arrays are compatible, but I can guarantee it. How can I do it?I have a tall array, that contains audio samples created from a dataset like this:
dataS = tall(fileDatastore(fullfile(appData.OutputDirectoryPath, ‘tmp_*.mat’), …
‘ReadFcn’, @(fname) getfield(load(fname), ‘resampled’), …
‘UniformRead’, true));
This data is large, but its element count is known. I need to multiply this data with a signal:
t = ((0:N-1)’/fsn);
carrier = exp(1j*2*pi*carrierFreq*t);
ccy = dataS.*carrier;
The last line fails, with the following error:
Error using .*
Incompatible non-scalar tall array arguments. Each of the tall arrays must be the same size in the first dimension, must be derived from a single tall array, and must not have been indexed differently in the first
dimension (indexing operations include functions such as VERTCAT, SPLITAPPLY, SORT, CELL2MAT, SYNCHRONIZE, RETIME and so on).
I assume the issue is that it cannot guarantee that the arrays are compatible, but I can guarantee it. How can I do it? I have a tall array, that contains audio samples created from a dataset like this:
dataS = tall(fileDatastore(fullfile(appData.OutputDirectoryPath, ‘tmp_*.mat’), …
‘ReadFcn’, @(fname) getfield(load(fname), ‘resampled’), …
‘UniformRead’, true));
This data is large, but its element count is known. I need to multiply this data with a signal:
t = ((0:N-1)’/fsn);
carrier = exp(1j*2*pi*carrierFreq*t);
ccy = dataS.*carrier;
The last line fails, with the following error:
Error using .*
Incompatible non-scalar tall array arguments. Each of the tall arrays must be the same size in the first dimension, must be derived from a single tall array, and must not have been indexed differently in the first
dimension (indexing operations include functions such as VERTCAT, SPLITAPPLY, SORT, CELL2MAT, SYNCHRONIZE, RETIME and so on).
I assume the issue is that it cannot guarantee that the arrays are compatible, but I can guarantee it. How can I do it? dataset, tall, signal processing MATLAB Answers — New Questions
How to interpolate and smooth across values in a matrix, while ignoring NaN values?
Hi there!
I am having some trouble figuring out how to interpolate across a matrix while ignoring NaN values. I have matrix data where each cell represents information from a single microphone, some of the microphones in the array malfunctioned, returning 0 values. I don’t want these microphones contributing 0s to my interpolation, as the 0s don’t represent real data, but rather nonfunctioning microphones. I set these malfunctioning microphones to NaN. My matrix is called ‘scan’ and looks something like this.
.3 .4 .3 .2 NaN .2
.4 .4 .3 .4 .2 .2
.3 NaN .5 .2 .4. .4
.3 .3 .4 .2 .2 .1
I need to interpolate across values in this matrix but I want the interpolation to ignore the NaN values, interpolating across them but without receiving them as input. Below is the interpolation code before I tried to have it ignore the NaN values
scan2 = scan’; %tranposes scan to scan2
interpImage = scatteredInterpolant(micsX(:),micsY(:), scan2(:), ‘natural’); %interpolates scan2 onto fine grid, calls it enInterp
enInterp = interpImage(meshImageX, meshImageY);
K = 0.125*ones(1, 1); %smooths the interpolated data
enInterpSmooth = conv2(enInterp, K, ‘valid’);
enInterpSmooth = enInterpSmooth’; %puts enInterpSmooth back to the original orientation
The above code returns the interpolation with huge chunks of NaNs surrounding where the NaN microphones were.
I added the below code to try and get an interpolation that ignored the NaN values, but it still returned large amounts of NaNs.
% Replace NaN values in enInterpSmooth with interpolated values
nanIndices = isnan(enInterpSmooth);
enInterpSmooth(nanIndices) = enInterp(nanIndices); % Replace NaNs with interpolated values
Any recommendations for how to remedy this problem would be much appreciated!!
Thank y’all so much!Hi there!
I am having some trouble figuring out how to interpolate across a matrix while ignoring NaN values. I have matrix data where each cell represents information from a single microphone, some of the microphones in the array malfunctioned, returning 0 values. I don’t want these microphones contributing 0s to my interpolation, as the 0s don’t represent real data, but rather nonfunctioning microphones. I set these malfunctioning microphones to NaN. My matrix is called ‘scan’ and looks something like this.
.3 .4 .3 .2 NaN .2
.4 .4 .3 .4 .2 .2
.3 NaN .5 .2 .4. .4
.3 .3 .4 .2 .2 .1
I need to interpolate across values in this matrix but I want the interpolation to ignore the NaN values, interpolating across them but without receiving them as input. Below is the interpolation code before I tried to have it ignore the NaN values
scan2 = scan’; %tranposes scan to scan2
interpImage = scatteredInterpolant(micsX(:),micsY(:), scan2(:), ‘natural’); %interpolates scan2 onto fine grid, calls it enInterp
enInterp = interpImage(meshImageX, meshImageY);
K = 0.125*ones(1, 1); %smooths the interpolated data
enInterpSmooth = conv2(enInterp, K, ‘valid’);
enInterpSmooth = enInterpSmooth’; %puts enInterpSmooth back to the original orientation
The above code returns the interpolation with huge chunks of NaNs surrounding where the NaN microphones were.
I added the below code to try and get an interpolation that ignored the NaN values, but it still returned large amounts of NaNs.
% Replace NaN values in enInterpSmooth with interpolated values
nanIndices = isnan(enInterpSmooth);
enInterpSmooth(nanIndices) = enInterp(nanIndices); % Replace NaNs with interpolated values
Any recommendations for how to remedy this problem would be much appreciated!!
Thank y’all so much! Hi there!
I am having some trouble figuring out how to interpolate across a matrix while ignoring NaN values. I have matrix data where each cell represents information from a single microphone, some of the microphones in the array malfunctioned, returning 0 values. I don’t want these microphones contributing 0s to my interpolation, as the 0s don’t represent real data, but rather nonfunctioning microphones. I set these malfunctioning microphones to NaN. My matrix is called ‘scan’ and looks something like this.
.3 .4 .3 .2 NaN .2
.4 .4 .3 .4 .2 .2
.3 NaN .5 .2 .4. .4
.3 .3 .4 .2 .2 .1
I need to interpolate across values in this matrix but I want the interpolation to ignore the NaN values, interpolating across them but without receiving them as input. Below is the interpolation code before I tried to have it ignore the NaN values
scan2 = scan’; %tranposes scan to scan2
interpImage = scatteredInterpolant(micsX(:),micsY(:), scan2(:), ‘natural’); %interpolates scan2 onto fine grid, calls it enInterp
enInterp = interpImage(meshImageX, meshImageY);
K = 0.125*ones(1, 1); %smooths the interpolated data
enInterpSmooth = conv2(enInterp, K, ‘valid’);
enInterpSmooth = enInterpSmooth’; %puts enInterpSmooth back to the original orientation
The above code returns the interpolation with huge chunks of NaNs surrounding where the NaN microphones were.
I added the below code to try and get an interpolation that ignored the NaN values, but it still returned large amounts of NaNs.
% Replace NaN values in enInterpSmooth with interpolated values
nanIndices = isnan(enInterpSmooth);
enInterpSmooth(nanIndices) = enInterp(nanIndices); % Replace NaNs with interpolated values
Any recommendations for how to remedy this problem would be much appreciated!!
Thank y’all so much! matrix, interpolation, nan, scatteredinterpolant, matrices MATLAB Answers — New Questions
Can I use the Open In MATLAB Online from MATLAB Grader feature across different browsers/machines at the same time?
I am sometimes loggied into multipe computers, or mutliple browsers on the same computer, using the same MathWorks account and running a MATLAB Online session in my browser. When I Open a MATLAB Grader problem in MATLAB Online, sometimes the MATLAB Grader problem content or left side panel either loads incompletely or doesn’t load at all. How can I get MATLAB Grader to launch successfully in MATLAB Online in this scenario?I am sometimes loggied into multipe computers, or mutliple browsers on the same computer, using the same MathWorks account and running a MATLAB Online session in my browser. When I Open a MATLAB Grader problem in MATLAB Online, sometimes the MATLAB Grader problem content or left side panel either loads incompletely or doesn’t load at all. How can I get MATLAB Grader to launch successfully in MATLAB Online in this scenario? I am sometimes loggied into multipe computers, or mutliple browsers on the same computer, using the same MathWorks account and running a MATLAB Online session in my browser. When I Open a MATLAB Grader problem in MATLAB Online, sometimes the MATLAB Grader problem content or left side panel either loads incompletely or doesn’t load at all. How can I get MATLAB Grader to launch successfully in MATLAB Online in this scenario? matlab_grader MATLAB Answers — New Questions
Set up a gradient descent algorithm for a multivariable economic dispatch problem
Hi,
I’m working on an electrical engineering optimization problem. I tried running the problem following an algorithm which was given to me in class, but I am getting some undefined results. I am attaching my code below:
clc; clear;
syms x y z
vX = [x y z].’; % set up x vector and transposition
Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]’; % calculate Q and b values according to quadratic problem formula, reproduced as f below
tau = 0.00001; % set up tolerance value
f = vX.’*Q*vX+b.’*vX+949;
Delf = [diff(f,x); diff(f,y); diff(f,z)];
t = (Delf.’*Delf)/(Delf.’*Q*Delf);
vx_old = [300 300 100]’; % set up initial vX value
grad = double(subs(Delf,vX,vx_old));
lenGrad = norm(grad);
Data = [vx_old;lenGrad];
while lenGrad > tau
t = subs(t,vX,vx_old);
vx_new = vx_old-t*grad;
grad = double(subs(Delf,vX,vx_new));
lenGrad = norm(grad);
Data = [Data [vx_new;lenGrad]];
vx_old = vx_new;
end
vx_new = double(vx_new)
fval = double(subs(f,vX,vx_new));
disp([‘Minimum value of the function is: ‘,num2str(fval)]);
disp([‘argmin of x value is: ‘,num2str(vx_new(1))]);
disp([‘argmin of y value is: ‘,num2str(vx_new(2))]);
disp([‘argmin of z value is: ‘,num2str(vx_new(3))]);
For reference, the objective function, which I am trying to minimize, is:
0.001562P1^2 + 0.00194P2^2 + 0.00482P3^2 + 7.92P1 + 7.85P2 + 7.97P3 + 949,
where each P represents the power output of 3 generating thermal units. I understand other methods are better suited to this problem, but I am required to try out unconstrained optimization on this problem. Would you help me determine if my code is wrong at some point or what I can do to improve it?Hi,
I’m working on an electrical engineering optimization problem. I tried running the problem following an algorithm which was given to me in class, but I am getting some undefined results. I am attaching my code below:
clc; clear;
syms x y z
vX = [x y z].’; % set up x vector and transposition
Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]’; % calculate Q and b values according to quadratic problem formula, reproduced as f below
tau = 0.00001; % set up tolerance value
f = vX.’*Q*vX+b.’*vX+949;
Delf = [diff(f,x); diff(f,y); diff(f,z)];
t = (Delf.’*Delf)/(Delf.’*Q*Delf);
vx_old = [300 300 100]’; % set up initial vX value
grad = double(subs(Delf,vX,vx_old));
lenGrad = norm(grad);
Data = [vx_old;lenGrad];
while lenGrad > tau
t = subs(t,vX,vx_old);
vx_new = vx_old-t*grad;
grad = double(subs(Delf,vX,vx_new));
lenGrad = norm(grad);
Data = [Data [vx_new;lenGrad]];
vx_old = vx_new;
end
vx_new = double(vx_new)
fval = double(subs(f,vX,vx_new));
disp([‘Minimum value of the function is: ‘,num2str(fval)]);
disp([‘argmin of x value is: ‘,num2str(vx_new(1))]);
disp([‘argmin of y value is: ‘,num2str(vx_new(2))]);
disp([‘argmin of z value is: ‘,num2str(vx_new(3))]);
For reference, the objective function, which I am trying to minimize, is:
0.001562P1^2 + 0.00194P2^2 + 0.00482P3^2 + 7.92P1 + 7.85P2 + 7.97P3 + 949,
where each P represents the power output of 3 generating thermal units. I understand other methods are better suited to this problem, but I am required to try out unconstrained optimization on this problem. Would you help me determine if my code is wrong at some point or what I can do to improve it? Hi,
I’m working on an electrical engineering optimization problem. I tried running the problem following an algorithm which was given to me in class, but I am getting some undefined results. I am attaching my code below:
clc; clear;
syms x y z
vX = [x y z].’; % set up x vector and transposition
Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]’; % calculate Q and b values according to quadratic problem formula, reproduced as f below
tau = 0.00001; % set up tolerance value
f = vX.’*Q*vX+b.’*vX+949;
Delf = [diff(f,x); diff(f,y); diff(f,z)];
t = (Delf.’*Delf)/(Delf.’*Q*Delf);
vx_old = [300 300 100]’; % set up initial vX value
grad = double(subs(Delf,vX,vx_old));
lenGrad = norm(grad);
Data = [vx_old;lenGrad];
while lenGrad > tau
t = subs(t,vX,vx_old);
vx_new = vx_old-t*grad;
grad = double(subs(Delf,vX,vx_new));
lenGrad = norm(grad);
Data = [Data [vx_new;lenGrad]];
vx_old = vx_new;
end
vx_new = double(vx_new)
fval = double(subs(f,vX,vx_new));
disp([‘Minimum value of the function is: ‘,num2str(fval)]);
disp([‘argmin of x value is: ‘,num2str(vx_new(1))]);
disp([‘argmin of y value is: ‘,num2str(vx_new(2))]);
disp([‘argmin of z value is: ‘,num2str(vx_new(3))]);
For reference, the objective function, which I am trying to minimize, is:
0.001562P1^2 + 0.00194P2^2 + 0.00482P3^2 + 7.92P1 + 7.85P2 + 7.97P3 + 949,
where each P represents the power output of 3 generating thermal units. I understand other methods are better suited to this problem, but I am required to try out unconstrained optimization on this problem. Would you help me determine if my code is wrong at some point or what I can do to improve it? optimization, gradient search, economic dispatch MATLAB Answers — New Questions
Export Simscape components to FMU
Hi,
I am currently trying to export single Simscape components as FMUs.
How do i create an interface to e.g. the physical water flow inside a pipe? I can only interface to Simulink signals, which can be converted to Simscape physical signals, but this doesn’t help me. I need to create an interface directly to the specific Simscape domain, e.g. water, current, air,….
If i create PMC ports on root level, no interface is created for those in the FMU. There is barely any documentation on how to export Simscape models to FMU, especially single components.
I am aware I can export a whole Simscape model that later works as a black box and only interfaces via Simulink signals, but this doesn’t help.
Best regards,
ChristianHi,
I am currently trying to export single Simscape components as FMUs.
How do i create an interface to e.g. the physical water flow inside a pipe? I can only interface to Simulink signals, which can be converted to Simscape physical signals, but this doesn’t help me. I need to create an interface directly to the specific Simscape domain, e.g. water, current, air,….
If i create PMC ports on root level, no interface is created for those in the FMU. There is barely any documentation on how to export Simscape models to FMU, especially single components.
I am aware I can export a whole Simscape model that later works as a black box and only interfaces via Simulink signals, but this doesn’t help.
Best regards,
Christian Hi,
I am currently trying to export single Simscape components as FMUs.
How do i create an interface to e.g. the physical water flow inside a pipe? I can only interface to Simulink signals, which can be converted to Simscape physical signals, but this doesn’t help me. I need to create an interface directly to the specific Simscape domain, e.g. water, current, air,….
If i create PMC ports on root level, no interface is created for those in the FMU. There is barely any documentation on how to export Simscape models to FMU, especially single components.
I am aware I can export a whole Simscape model that later works as a black box and only interfaces via Simulink signals, but this doesn’t help.
Best regards,
Christian simscape fmu export MATLAB Answers — New Questions