Author: PuTI
What is the correct way to save a large MATLAB structure?
I have a MATLAB structure which is just over 21GB in memory (from whos) and when I save this to a MAT file with the "-v7.3" and "-nocompression" flags it takes well over an hour (on a high performance workstation with an NVMe SSD) and I get a file which is 77GB on disk. I understand that there is some overhead in saving to a MAT file and that "-nocompression" will result in larger files that with compression (but I gave up after about 3 hours waiting for a compressed version to save), but how can 56GB of "overhead" be considered acceptable?
I only need to save this one structure and I won’t be adding any other data or modifying the MAT file, so all of the additional features of the v7.3 format are of no use to me, I just need support for >2Gb variables. I attempted to use the undocumented getByteStreamFromArray function to get a byte array I can just dump to a file but this just returned "Error during serialization".
Am I somehow missing a "correct" way to do this efficiently? Or are my only options to either split my data in to a bunch of <2Gb variables to save in a v7 format or write my own serializer? I appreciate that doing either of these isn’t exactly a massive job, I’m just very surprised there isn’t better native support for large files!I have a MATLAB structure which is just over 21GB in memory (from whos) and when I save this to a MAT file with the "-v7.3" and "-nocompression" flags it takes well over an hour (on a high performance workstation with an NVMe SSD) and I get a file which is 77GB on disk. I understand that there is some overhead in saving to a MAT file and that "-nocompression" will result in larger files that with compression (but I gave up after about 3 hours waiting for a compressed version to save), but how can 56GB of "overhead" be considered acceptable?
I only need to save this one structure and I won’t be adding any other data or modifying the MAT file, so all of the additional features of the v7.3 format are of no use to me, I just need support for >2Gb variables. I attempted to use the undocumented getByteStreamFromArray function to get a byte array I can just dump to a file but this just returned "Error during serialization".
Am I somehow missing a "correct" way to do this efficiently? Or are my only options to either split my data in to a bunch of <2Gb variables to save in a v7 format or write my own serializer? I appreciate that doing either of these isn’t exactly a massive job, I’m just very surprised there isn’t better native support for large files! I have a MATLAB structure which is just over 21GB in memory (from whos) and when I save this to a MAT file with the "-v7.3" and "-nocompression" flags it takes well over an hour (on a high performance workstation with an NVMe SSD) and I get a file which is 77GB on disk. I understand that there is some overhead in saving to a MAT file and that "-nocompression" will result in larger files that with compression (but I gave up after about 3 hours waiting for a compressed version to save), but how can 56GB of "overhead" be considered acceptable?
I only need to save this one structure and I won’t be adding any other data or modifying the MAT file, so all of the additional features of the v7.3 format are of no use to me, I just need support for >2Gb variables. I attempted to use the undocumented getByteStreamFromArray function to get a byte array I can just dump to a file but this just returned "Error during serialization".
Am I somehow missing a "correct" way to do this efficiently? Or are my only options to either split my data in to a bunch of <2Gb variables to save in a v7 format or write my own serializer? I appreciate that doing either of these isn’t exactly a massive job, I’m just very surprised there isn’t better native support for large files! save, big-data MATLAB Answers — New Questions
Create a Cartesian product with user-specified constraints
Hello,
I’m trying to figure out how to code a simple and flexible cartesian product function
To abstract how it works conceptually, let’s imagine we have a 4 X 10 matrix, with each column having ordered values 1 through 4, like this
1 1 1 1 1 1 etc
2 2 2 2 2 2 etc
3 3 3 3 3 3 etc
4 4 4 4 4 4 etc
The user specifies the number of rows and columns to be used to make the cartesian product. For example, if they specify 2 rows and 3 columns, the result should be (although the order doesn’t matter)
1 1 1
1 1 2
1 2 1
2 1 1
1 2 2
2 1 2
2 2 1
2 2 2
The closest function that I’ve found is
function C = cartesian(varargin)
% This function creates cartesian products from input arrays
args = varargin;
n = nargin;
[F{1:n}] = ndgrid(args{:});
for i=n:-1:1
G(:,i) = F{i}(:);
end
C = unique(G , ‘rows’);
end
However, it asks for a series of arrays as input. In my case, I need the input to be a matrix of variable dimensions. Thank you for any help!Hello,
I’m trying to figure out how to code a simple and flexible cartesian product function
To abstract how it works conceptually, let’s imagine we have a 4 X 10 matrix, with each column having ordered values 1 through 4, like this
1 1 1 1 1 1 etc
2 2 2 2 2 2 etc
3 3 3 3 3 3 etc
4 4 4 4 4 4 etc
The user specifies the number of rows and columns to be used to make the cartesian product. For example, if they specify 2 rows and 3 columns, the result should be (although the order doesn’t matter)
1 1 1
1 1 2
1 2 1
2 1 1
1 2 2
2 1 2
2 2 1
2 2 2
The closest function that I’ve found is
function C = cartesian(varargin)
% This function creates cartesian products from input arrays
args = varargin;
n = nargin;
[F{1:n}] = ndgrid(args{:});
for i=n:-1:1
G(:,i) = F{i}(:);
end
C = unique(G , ‘rows’);
end
However, it asks for a series of arrays as input. In my case, I need the input to be a matrix of variable dimensions. Thank you for any help! Hello,
I’m trying to figure out how to code a simple and flexible cartesian product function
To abstract how it works conceptually, let’s imagine we have a 4 X 10 matrix, with each column having ordered values 1 through 4, like this
1 1 1 1 1 1 etc
2 2 2 2 2 2 etc
3 3 3 3 3 3 etc
4 4 4 4 4 4 etc
The user specifies the number of rows and columns to be used to make the cartesian product. For example, if they specify 2 rows and 3 columns, the result should be (although the order doesn’t matter)
1 1 1
1 1 2
1 2 1
2 1 1
1 2 2
2 1 2
2 2 1
2 2 2
The closest function that I’ve found is
function C = cartesian(varargin)
% This function creates cartesian products from input arrays
args = varargin;
n = nargin;
[F{1:n}] = ndgrid(args{:});
for i=n:-1:1
G(:,i) = F{i}(:);
end
C = unique(G , ‘rows’);
end
However, it asks for a series of arrays as input. In my case, I need the input to be a matrix of variable dimensions. Thank you for any help! cartesian product MATLAB Answers — New Questions
Find model of a 2 to 4000 horsepower squirrel cage induction motor on a variable speed drive
Speed regulators are a poor choice for Mine Hoist control. Mine hoists consist of distributed masses connected with very long springy ropes tied to the conveyances. Thus the system can present varying natural frequencies. ALSO Double drum hoists present a mass to the speed regulator that varies by a factor or 4 or 5 times. Hoist clutched in (two drums full of rope) to clutched out (single drum no rope)
A much better choice: The speed torque characteristics of large synchronous and induction motors are well suited for this application. Mine Hoists are stand alone machines that do not require the speed and accuracy required in steel and paper mills.
I would like to validate these statments using Matlab simulation.Speed regulators are a poor choice for Mine Hoist control. Mine hoists consist of distributed masses connected with very long springy ropes tied to the conveyances. Thus the system can present varying natural frequencies. ALSO Double drum hoists present a mass to the speed regulator that varies by a factor or 4 or 5 times. Hoist clutched in (two drums full of rope) to clutched out (single drum no rope)
A much better choice: The speed torque characteristics of large synchronous and induction motors are well suited for this application. Mine Hoists are stand alone machines that do not require the speed and accuracy required in steel and paper mills.
I would like to validate these statments using Matlab simulation. Speed regulators are a poor choice for Mine Hoist control. Mine hoists consist of distributed masses connected with very long springy ropes tied to the conveyances. Thus the system can present varying natural frequencies. ALSO Double drum hoists present a mass to the speed regulator that varies by a factor or 4 or 5 times. Hoist clutched in (two drums full of rope) to clutched out (single drum no rope)
A much better choice: The speed torque characteristics of large synchronous and induction motors are well suited for this application. Mine Hoists are stand alone machines that do not require the speed and accuracy required in steel and paper mills.
I would like to validate these statments using Matlab simulation. asynchrous motor, sycronous motor, motor characteristics, variable speed drive, mine hoist MATLAB Answers — New Questions
How to use the overpotential in the parameterised battery for the simscape battery equivalent circuit? I am trying to use the molicell 21700 cell. I can run the R0 model.
The simscape equivalent circuit battery has a bunch of parameterized cells. I am using the molicelI NR_21700_P45B lithium ion battery. The in built data set has 101 state of charge breakpoints, and 7 temperature levels. So, the R0 value has a matrix of size 101×7. However, the overpotential (R1-C1, R2-C2, R3-C3) has the size 7×3. Obviously, when I run the model, I get the following error :
size of First polarization resistance, R1(SOC,T) must be equal to length of State of charge breakpoints for resistance, SOC by length of Temperature breakpoints for resistance, T.
I have read the model documentation clearly, there is a part where a hint of how to use the values properly is given: The software does not parameterize the change in series resistance with the battery cycle life and the dynamic RC network parameters. Instead, the software sums the net resistance of the RC network resistors to the series resistance of the specific parameterized data. To populate the RC parameter data, subtract the net RC network resistance from the series resistance data.
This description to me is not very clear. The model contains the following data heads: Instantaneous resistance, R0 (SOC, T), first polarization resistance R1(SOC, T) and Tau, if I use one time-constant dynamics, which are of my interest. There is nothing mentioned explicitly about a series resistance from which the net RC resistance is to be subtracted. A demonstration of how to use this will be useful.The simscape equivalent circuit battery has a bunch of parameterized cells. I am using the molicelI NR_21700_P45B lithium ion battery. The in built data set has 101 state of charge breakpoints, and 7 temperature levels. So, the R0 value has a matrix of size 101×7. However, the overpotential (R1-C1, R2-C2, R3-C3) has the size 7×3. Obviously, when I run the model, I get the following error :
size of First polarization resistance, R1(SOC,T) must be equal to length of State of charge breakpoints for resistance, SOC by length of Temperature breakpoints for resistance, T.
I have read the model documentation clearly, there is a part where a hint of how to use the values properly is given: The software does not parameterize the change in series resistance with the battery cycle life and the dynamic RC network parameters. Instead, the software sums the net resistance of the RC network resistors to the series resistance of the specific parameterized data. To populate the RC parameter data, subtract the net RC network resistance from the series resistance data.
This description to me is not very clear. The model contains the following data heads: Instantaneous resistance, R0 (SOC, T), first polarization resistance R1(SOC, T) and Tau, if I use one time-constant dynamics, which are of my interest. There is nothing mentioned explicitly about a series resistance from which the net RC resistance is to be subtracted. A demonstration of how to use this will be useful. The simscape equivalent circuit battery has a bunch of parameterized cells. I am using the molicelI NR_21700_P45B lithium ion battery. The in built data set has 101 state of charge breakpoints, and 7 temperature levels. So, the R0 value has a matrix of size 101×7. However, the overpotential (R1-C1, R2-C2, R3-C3) has the size 7×3. Obviously, when I run the model, I get the following error :
size of First polarization resistance, R1(SOC,T) must be equal to length of State of charge breakpoints for resistance, SOC by length of Temperature breakpoints for resistance, T.
I have read the model documentation clearly, there is a part where a hint of how to use the values properly is given: The software does not parameterize the change in series resistance with the battery cycle life and the dynamic RC network parameters. Instead, the software sums the net resistance of the RC network resistors to the series resistance of the specific parameterized data. To populate the RC parameter data, subtract the net RC network resistance from the series resistance data.
This description to me is not very clear. The model contains the following data heads: Instantaneous resistance, R0 (SOC, T), first polarization resistance R1(SOC, T) and Tau, if I use one time-constant dynamics, which are of my interest. There is nothing mentioned explicitly about a series resistance from which the net RC resistance is to be subtracted. A demonstration of how to use this will be useful. simulink, simscape, matrices, battery, equivalent_circuit, lithium-ion, cell modelling, matrix manipulation MATLAB Answers — New Questions
[Vehicle Dynamics Blockset] Confused about coordinate system conventions
Hello everyone!
Recently I’m working on this Simulink example: https://www.mathworks.com/help/vdynblks/ug/braking-test.html
I also checked the setting of coordinate system in VDBS: https://www.mathworks.com/help/vdynblks/ug/coordinate-systems-in-vehicle-dynamics-blockset.html
In this simulation, I assume a few settings:
(1) the ground is completely flat and horizontal
(2) the origin of inertial (world-fixed) frame is attach at the ground surface, with Z-axe pointing downward.
(3) the origin of tire frame is attach at the ground surface, with Z-axe pointing upward.
I want to slightly modify the Simulink example shared above, by changing parameters in some blocks:
(1) Block 1 —- Vehicle Body 6DOF
Based on the figure below, the vehicle-fixed frame should attach at the center of mass (CoM)
In this block, I want to change the initial position of CoM in inertial frame (aka, world-fixed frame), as shown below:
Since the inertial frame is attached on the ground, for the initial vertical position of CoM, I set it to -0.30938-VEH.HeightCG = -0.44338. Here, the value 0.30938 is the unload tire radius, and VEH.HeightCG is the vertical offset between CoM and axle plane (h in figure "vehicle layout" above). From my understanding, this position can let the tire just slightly contact the ground.
(2) Block 2 —- Combined Slip Wheel 2DOF
In this block, I assign a small value to this parameter:
Based on my understanding, this value refers to the compression of tire due to vertical load. A positive value means that the tire is compressed.
(3) Block 3 —- Ground Feedback
Since I assume inertial frame is attach at the ground surface, so the ground height for all wheels G_xx_z (xx: FL, FR, RL, RR) are set to 0 all the time.
Then I start the simulation. During the first 2 seconds, the vehicle will stablize on the ground. I logged some relavent signals to check the vehicle states:
(1) Passenger Vehicle:1.Body.InertFrm.Cg.Disp.Z —- the vertical position of CoM in inertial frame
(2) Passenger Vehicle:1.Body.InertFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in inertial frame
(3) Passenger Vehicle:1.Body.BdyFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in vehicle-fixed frame
(4) Passenger Vehicle:1.Wheels.TireFrame.z(1) —- the vertical displacement of front-left axle in tire frame? (I’m not sure, just based on the description in https://www.mathworks.com/help/vdynblks/ref/combinedslipwheel2dof.html)
(5) Passenger Vehicle:1.Wheels.TireFrame.dz(1) —- the vertical displacement of front-left axle (wheel) in vehicle-fixed frame? (I’m not sure, just based on the computation in the model, as shown below)
The stablizing process for vehicle is shown below:
After the vehicle is stablized on the ground, we can see:
(1) the vertical position of CoM in inertial frame converges to around 0.0
(2) the vertical position of front-left axle in inertial frame is close to that value in vehicle-fixed frame, while the latter one is always a constant value equal to VEH.HeightCG
(3) TireFrame.z(1) and TireFrame.dz(1) converge to some small values near 0.0
(4) The initial value of TireFrame.z(1) is the opposite of the initial value of a parameter z0 in Block "Combined Slip Wheel 2DOF" (as shown above). It might be due to a flip of Z-axe direction.
It means that the CoM of the vehicle is near the ground, which means that the vehicle chassis penetrates and sinks under the ground. It doesn’t make sense to me. I have a few questions:
(1) What are the setting of coordinate systems for different components (CoM, axles, wheels, tires)? —- where should the origin of a specific frame attach?
(2) What are the exact meaning of TireFrame.z and TireFrame.dz? The vertical position of wheel center? In which frame?
I would be very grateful if someone could provide an answer.
Regards,
Zijun HeHello everyone!
Recently I’m working on this Simulink example: https://www.mathworks.com/help/vdynblks/ug/braking-test.html
I also checked the setting of coordinate system in VDBS: https://www.mathworks.com/help/vdynblks/ug/coordinate-systems-in-vehicle-dynamics-blockset.html
In this simulation, I assume a few settings:
(1) the ground is completely flat and horizontal
(2) the origin of inertial (world-fixed) frame is attach at the ground surface, with Z-axe pointing downward.
(3) the origin of tire frame is attach at the ground surface, with Z-axe pointing upward.
I want to slightly modify the Simulink example shared above, by changing parameters in some blocks:
(1) Block 1 —- Vehicle Body 6DOF
Based on the figure below, the vehicle-fixed frame should attach at the center of mass (CoM)
In this block, I want to change the initial position of CoM in inertial frame (aka, world-fixed frame), as shown below:
Since the inertial frame is attached on the ground, for the initial vertical position of CoM, I set it to -0.30938-VEH.HeightCG = -0.44338. Here, the value 0.30938 is the unload tire radius, and VEH.HeightCG is the vertical offset between CoM and axle plane (h in figure "vehicle layout" above). From my understanding, this position can let the tire just slightly contact the ground.
(2) Block 2 —- Combined Slip Wheel 2DOF
In this block, I assign a small value to this parameter:
Based on my understanding, this value refers to the compression of tire due to vertical load. A positive value means that the tire is compressed.
(3) Block 3 —- Ground Feedback
Since I assume inertial frame is attach at the ground surface, so the ground height for all wheels G_xx_z (xx: FL, FR, RL, RR) are set to 0 all the time.
Then I start the simulation. During the first 2 seconds, the vehicle will stablize on the ground. I logged some relavent signals to check the vehicle states:
(1) Passenger Vehicle:1.Body.InertFrm.Cg.Disp.Z —- the vertical position of CoM in inertial frame
(2) Passenger Vehicle:1.Body.InertFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in inertial frame
(3) Passenger Vehicle:1.Body.BdyFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in vehicle-fixed frame
(4) Passenger Vehicle:1.Wheels.TireFrame.z(1) —- the vertical displacement of front-left axle in tire frame? (I’m not sure, just based on the description in https://www.mathworks.com/help/vdynblks/ref/combinedslipwheel2dof.html)
(5) Passenger Vehicle:1.Wheels.TireFrame.dz(1) —- the vertical displacement of front-left axle (wheel) in vehicle-fixed frame? (I’m not sure, just based on the computation in the model, as shown below)
The stablizing process for vehicle is shown below:
After the vehicle is stablized on the ground, we can see:
(1) the vertical position of CoM in inertial frame converges to around 0.0
(2) the vertical position of front-left axle in inertial frame is close to that value in vehicle-fixed frame, while the latter one is always a constant value equal to VEH.HeightCG
(3) TireFrame.z(1) and TireFrame.dz(1) converge to some small values near 0.0
(4) The initial value of TireFrame.z(1) is the opposite of the initial value of a parameter z0 in Block "Combined Slip Wheel 2DOF" (as shown above). It might be due to a flip of Z-axe direction.
It means that the CoM of the vehicle is near the ground, which means that the vehicle chassis penetrates and sinks under the ground. It doesn’t make sense to me. I have a few questions:
(1) What are the setting of coordinate systems for different components (CoM, axles, wheels, tires)? —- where should the origin of a specific frame attach?
(2) What are the exact meaning of TireFrame.z and TireFrame.dz? The vertical position of wheel center? In which frame?
I would be very grateful if someone could provide an answer.
Regards,
Zijun He Hello everyone!
Recently I’m working on this Simulink example: https://www.mathworks.com/help/vdynblks/ug/braking-test.html
I also checked the setting of coordinate system in VDBS: https://www.mathworks.com/help/vdynblks/ug/coordinate-systems-in-vehicle-dynamics-blockset.html
In this simulation, I assume a few settings:
(1) the ground is completely flat and horizontal
(2) the origin of inertial (world-fixed) frame is attach at the ground surface, with Z-axe pointing downward.
(3) the origin of tire frame is attach at the ground surface, with Z-axe pointing upward.
I want to slightly modify the Simulink example shared above, by changing parameters in some blocks:
(1) Block 1 —- Vehicle Body 6DOF
Based on the figure below, the vehicle-fixed frame should attach at the center of mass (CoM)
In this block, I want to change the initial position of CoM in inertial frame (aka, world-fixed frame), as shown below:
Since the inertial frame is attached on the ground, for the initial vertical position of CoM, I set it to -0.30938-VEH.HeightCG = -0.44338. Here, the value 0.30938 is the unload tire radius, and VEH.HeightCG is the vertical offset between CoM and axle plane (h in figure "vehicle layout" above). From my understanding, this position can let the tire just slightly contact the ground.
(2) Block 2 —- Combined Slip Wheel 2DOF
In this block, I assign a small value to this parameter:
Based on my understanding, this value refers to the compression of tire due to vertical load. A positive value means that the tire is compressed.
(3) Block 3 —- Ground Feedback
Since I assume inertial frame is attach at the ground surface, so the ground height for all wheels G_xx_z (xx: FL, FR, RL, RR) are set to 0 all the time.
Then I start the simulation. During the first 2 seconds, the vehicle will stablize on the ground. I logged some relavent signals to check the vehicle states:
(1) Passenger Vehicle:1.Body.InertFrm.Cg.Disp.Z —- the vertical position of CoM in inertial frame
(2) Passenger Vehicle:1.Body.InertFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in inertial frame
(3) Passenger Vehicle:1.Body.BdyFrm.FrntAxl.Lft.Disp.Z —- the vertical position of front-left axle in vehicle-fixed frame
(4) Passenger Vehicle:1.Wheels.TireFrame.z(1) —- the vertical displacement of front-left axle in tire frame? (I’m not sure, just based on the description in https://www.mathworks.com/help/vdynblks/ref/combinedslipwheel2dof.html)
(5) Passenger Vehicle:1.Wheels.TireFrame.dz(1) —- the vertical displacement of front-left axle (wheel) in vehicle-fixed frame? (I’m not sure, just based on the computation in the model, as shown below)
The stablizing process for vehicle is shown below:
After the vehicle is stablized on the ground, we can see:
(1) the vertical position of CoM in inertial frame converges to around 0.0
(2) the vertical position of front-left axle in inertial frame is close to that value in vehicle-fixed frame, while the latter one is always a constant value equal to VEH.HeightCG
(3) TireFrame.z(1) and TireFrame.dz(1) converge to some small values near 0.0
(4) The initial value of TireFrame.z(1) is the opposite of the initial value of a parameter z0 in Block "Combined Slip Wheel 2DOF" (as shown above). It might be due to a flip of Z-axe direction.
It means that the CoM of the vehicle is near the ground, which means that the vehicle chassis penetrates and sinks under the ground. It doesn’t make sense to me. I have a few questions:
(1) What are the setting of coordinate systems for different components (CoM, axles, wheels, tires)? —- where should the origin of a specific frame attach?
(2) What are the exact meaning of TireFrame.z and TireFrame.dz? The vertical position of wheel center? In which frame?
I would be very grateful if someone could provide an answer.
Regards,
Zijun He vehicle dynamcis, coordinate system, vehicle dynamics blockset MATLAB Answers — New Questions
R2025b Update 4 does not support element-wise division (./) for tf objects?
The following should divide 1/(s+1) by 2, and 1/(s+2) by 3:
s = tf(‘s’);
sys1 = [1/(s+1), 1/(s+2)];
sys2 = [2, 3];
result = sys1 ./ sys2;
However, I get the following error message:
Error using ./
Invalid data type. Argument must be numeric, char, or logical.
Error in
test (line 4)
result = sys1 ./ sys2;The following should divide 1/(s+1) by 2, and 1/(s+2) by 3:
s = tf(‘s’);
sys1 = [1/(s+1), 1/(s+2)];
sys2 = [2, 3];
result = sys1 ./ sys2;
However, I get the following error message:
Error using ./
Invalid data type. Argument must be numeric, char, or logical.
Error in
test (line 4)
result = sys1 ./ sys2; The following should divide 1/(s+1) by 2, and 1/(s+2) by 3:
s = tf(‘s’);
sys1 = [1/(s+1), 1/(s+2)];
sys2 = [2, 3];
result = sys1 ./ sys2;
However, I get the following error message:
Error using ./
Invalid data type. Argument must be numeric, char, or logical.
Error in
test (line 4)
result = sys1 ./ sys2; control toolbox, element-wise division, transfer-function object MATLAB Answers — New Questions
How to draw the axis in the origin of the coordinate system?
Hello,
I want to make a figure that looks like this:
The code has always worked fine, but now it does not save the figure correctly anymore. (Maybe I updated MATLAB since last year, but I don’t know; I am using R2025b and don’t know, which one I used before.)
The exported figure looks like this:
The labels for f(x) and x are at the wrong places.
The error happens with exportfigure, print and manually saving the figure and png, jpg and pdf.
ChatGPT tells me, the problem is with "origin" but that there are only manual fixes which seem very complicated and unsatisfying to me.
Is there a fix for my figure / what should I do?
Thank you very much in advance!
%% Minimal Working Example
figureFontSize = 22;
set(0,’DefaultAxesFontSize’,22)
lwidth=2;
xMin=-5;
xMax=+6;
xVal=(xMin:0.05:xMax);
% Aufgabe
yVal(:,1)=-2*(xVal-4).^2+64;
% Aufgabe
fig=figure(1);
plot(xVal,yVal(:,1),’-k’,’LineWidth’,lwidth)
xlabel(‘x’)
ylA=ylabel(‘f(x)’);
xlim([xMin xMax])
a=max(abs(min(yVal(:,1))),max(yVal(:,1)));
ylim([-a-0.01*a +a])
set(gca, ‘xtick’, [])
set(gca, ‘ytick’, [])
ax = gca;
box(‘off’)
grid(‘off’)
set(gca, ‘xticklabel’, [])
set(gca, ‘yticklabel’, [])
ax.XAxisLocation = "origin";
ax.YAxisLocation = "origin";
% print(gcf,’test.png’,’-dpng’,’-r300′)
exportgraphics(gcf,’test.png’,’Resolution’,300)Hello,
I want to make a figure that looks like this:
The code has always worked fine, but now it does not save the figure correctly anymore. (Maybe I updated MATLAB since last year, but I don’t know; I am using R2025b and don’t know, which one I used before.)
The exported figure looks like this:
The labels for f(x) and x are at the wrong places.
The error happens with exportfigure, print and manually saving the figure and png, jpg and pdf.
ChatGPT tells me, the problem is with "origin" but that there are only manual fixes which seem very complicated and unsatisfying to me.
Is there a fix for my figure / what should I do?
Thank you very much in advance!
%% Minimal Working Example
figureFontSize = 22;
set(0,’DefaultAxesFontSize’,22)
lwidth=2;
xMin=-5;
xMax=+6;
xVal=(xMin:0.05:xMax);
% Aufgabe
yVal(:,1)=-2*(xVal-4).^2+64;
% Aufgabe
fig=figure(1);
plot(xVal,yVal(:,1),’-k’,’LineWidth’,lwidth)
xlabel(‘x’)
ylA=ylabel(‘f(x)’);
xlim([xMin xMax])
a=max(abs(min(yVal(:,1))),max(yVal(:,1)));
ylim([-a-0.01*a +a])
set(gca, ‘xtick’, [])
set(gca, ‘ytick’, [])
ax = gca;
box(‘off’)
grid(‘off’)
set(gca, ‘xticklabel’, [])
set(gca, ‘yticklabel’, [])
ax.XAxisLocation = "origin";
ax.YAxisLocation = "origin";
% print(gcf,’test.png’,’-dpng’,’-r300′)
exportgraphics(gcf,’test.png’,’Resolution’,300) Hello,
I want to make a figure that looks like this:
The code has always worked fine, but now it does not save the figure correctly anymore. (Maybe I updated MATLAB since last year, but I don’t know; I am using R2025b and don’t know, which one I used before.)
The exported figure looks like this:
The labels for f(x) and x are at the wrong places.
The error happens with exportfigure, print and manually saving the figure and png, jpg and pdf.
ChatGPT tells me, the problem is with "origin" but that there are only manual fixes which seem very complicated and unsatisfying to me.
Is there a fix for my figure / what should I do?
Thank you very much in advance!
%% Minimal Working Example
figureFontSize = 22;
set(0,’DefaultAxesFontSize’,22)
lwidth=2;
xMin=-5;
xMax=+6;
xVal=(xMin:0.05:xMax);
% Aufgabe
yVal(:,1)=-2*(xVal-4).^2+64;
% Aufgabe
fig=figure(1);
plot(xVal,yVal(:,1),’-k’,’LineWidth’,lwidth)
xlabel(‘x’)
ylA=ylabel(‘f(x)’);
xlim([xMin xMax])
a=max(abs(min(yVal(:,1))),max(yVal(:,1)));
ylim([-a-0.01*a +a])
set(gca, ‘xtick’, [])
set(gca, ‘ytick’, [])
ax = gca;
box(‘off’)
grid(‘off’)
set(gca, ‘xticklabel’, [])
set(gca, ‘yticklabel’, [])
ax.XAxisLocation = "origin";
ax.YAxisLocation = "origin";
% print(gcf,’test.png’,’-dpng’,’-r300′)
exportgraphics(gcf,’test.png’,’Resolution’,300) axislocation, origin MATLAB Answers — New Questions
How do I find if a point is within the volume of a rotated ellipsoid
I have an ellipsoid
[ex,ey,ez]=ellipsoid(long2,lat2,h2,bb2,aa2,cc2);
direction=[1 0 0];
s=surf(ex,ey,ez,’FaceAlpha’,0.05)
rotate(s,direction,theta(i),[long2 lat2 h2])
and say I have a point ;
point = [longX latX hX];
How do I find if point is within the volume of the rotated ellipsoid?
Thanks!I have an ellipsoid
[ex,ey,ez]=ellipsoid(long2,lat2,h2,bb2,aa2,cc2);
direction=[1 0 0];
s=surf(ex,ey,ez,’FaceAlpha’,0.05)
rotate(s,direction,theta(i),[long2 lat2 h2])
and say I have a point ;
point = [longX latX hX];
How do I find if point is within the volume of the rotated ellipsoid?
Thanks! I have an ellipsoid
[ex,ey,ez]=ellipsoid(long2,lat2,h2,bb2,aa2,cc2);
direction=[1 0 0];
s=surf(ex,ey,ez,’FaceAlpha’,0.05)
rotate(s,direction,theta(i),[long2 lat2 h2])
and say I have a point ;
point = [longX latX hX];
How do I find if point is within the volume of the rotated ellipsoid?
Thanks! ellipsoid, volume MATLAB Answers — New Questions
code improvement/optimization
I would like to know if this code can be improved or it well written? any comments will be greatly appreciated. Thanks
u0 = [1; 0; 0];
tspan = [0, 1e-3]; % small final time for explicit solver
tol = 1e-2;
h0 = 1e-6;
% Execute solver
[t, u, err, h_steps, stats] = dormand_prince_solver(@robertson_ode, tspan, u0, tol, h0);
% Part (a): Plot Solutions
figure;
subplot(2,1,1);
plot(t, u(:,1), ‘b’, t, u(:,3), ‘g’, ‘LineWidth’, 1.5);
legend(‘y1 (Reactant)’, ‘y3 (Product)’);
xlabel(‘Time t’); ylabel(‘Concentration’);
%title(‘Robertson Problem: Concentrations y1 and y3’);
grid on;
subplot(2,1,2);
semilogx(t, u(:,2), ‘r’, ‘LineWidth’, 1.5);
xlabel(‘Time t (Log Scale)’); ylabel(‘Concentration y2’);
%title(‘Robertson Problem: Intermediate Species y2’);
grid on;
% Part (b): Error and Stepsize Plots
figure;
subplot(2,1,1);
semilogy(t(2:end), err);
xlabel(‘Time t’); ylabel(‘||y_n^5 – y_n^4||_{infty} / h’);
%title(‘Part (b): Estimated Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps);
xlabel(‘Time t’); ylabel(‘Stepsize h’);
%title(‘Part (b): Stepsize h vs Time’);
grid on;
figure;
subplot(2,1,1);
semilogy(t(2:end), err, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Error’);
title(‘Estimated Local Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Stepsize h’);
title(‘Stepsize vs Time’);
grid on;
function dydt = robertson_ode(~, y)
% Coefficients from Screenshot 2026-02-28 at 3.38.10 PM.png
alpha = 0.04;
beta = 1e4;
gamma = 3e7;
dydt = [ -alpha*y(1) + beta*y(2)*y(3);
alpha*y(1) – beta*y(2)*y(3) – gamma*y(2)^2;
gamma*y(2)^2 ];
end
function [t_all, u_all, error_history, h_history, stats] = …
dormand_prince_solver(f, tspan, u0, tol, h0)
% Initialization
t = tspan(1);
u = u0(:);
h = h0;
t_all = t;
u_all = u.’;
error_history = [];
h_history = [];
n_accepted = 0;
n_rejected = 0;
% —- Dormand–Prince 4(5) Coefficients —-
a = [1/5, 0, 0, 0, 0, 0;
3/40, 9/40, 0, 0, 0, 0;
44/45, -56/15, 32/9, 0, 0, 0;
19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0;
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0;
35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
b4 = [5179/57600, 0, 7571/16695, 393/640, …
-92097/339200, 187/2100, 1/40];
b5 = [35/384, 0, 500/1113, 125/192, …
-2187/6784, 11/84, 0];
c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
safety = 0.9;
h_min = 1e-12;
h_max = 1.0;
while t < tspan(2)
if h < h_min
error(‘Step size underflow: problem likely stiff.’);
end
if t + h > tspan(2)
h = tspan(2) – t;
end
% —- Stage calculations —-
k = zeros(length(u), 7);
k(:,1) = f(t, u);
for i = 2:7
k(:,i) = f(t + c(i)*h, …
u + h * (k(:,1:i-1) * a(i-1,1:i-1)’));
end
% —- 4th and 5th order solutions —-
y4 = u + h * (k * b4′);
y5 = u + h * (k * b5′);
error_val = norm(y5 – y4, inf);
if error_val <= tol
t = t + h;
u = y5;
t_all = [t_all; t];
u_all = [u_all; u.’];
error_history = [error_history; error_val];
h_history = [h_history; h];
n_accepted = n_accepted + 1;
else
n_rejected = n_rejected + 1;
end
if error_val == 0
h_new = 2*h;
else
h_new = h * safety * (tol / error_val)^(1/5);
end
h = min(max(h_new, h_min), h_max);
end
stats.accepted = n_accepted;
stats.rejected = n_rejected;
endI would like to know if this code can be improved or it well written? any comments will be greatly appreciated. Thanks
u0 = [1; 0; 0];
tspan = [0, 1e-3]; % small final time for explicit solver
tol = 1e-2;
h0 = 1e-6;
% Execute solver
[t, u, err, h_steps, stats] = dormand_prince_solver(@robertson_ode, tspan, u0, tol, h0);
% Part (a): Plot Solutions
figure;
subplot(2,1,1);
plot(t, u(:,1), ‘b’, t, u(:,3), ‘g’, ‘LineWidth’, 1.5);
legend(‘y1 (Reactant)’, ‘y3 (Product)’);
xlabel(‘Time t’); ylabel(‘Concentration’);
%title(‘Robertson Problem: Concentrations y1 and y3’);
grid on;
subplot(2,1,2);
semilogx(t, u(:,2), ‘r’, ‘LineWidth’, 1.5);
xlabel(‘Time t (Log Scale)’); ylabel(‘Concentration y2’);
%title(‘Robertson Problem: Intermediate Species y2’);
grid on;
% Part (b): Error and Stepsize Plots
figure;
subplot(2,1,1);
semilogy(t(2:end), err);
xlabel(‘Time t’); ylabel(‘||y_n^5 – y_n^4||_{infty} / h’);
%title(‘Part (b): Estimated Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps);
xlabel(‘Time t’); ylabel(‘Stepsize h’);
%title(‘Part (b): Stepsize h vs Time’);
grid on;
figure;
subplot(2,1,1);
semilogy(t(2:end), err, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Error’);
title(‘Estimated Local Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Stepsize h’);
title(‘Stepsize vs Time’);
grid on;
function dydt = robertson_ode(~, y)
% Coefficients from Screenshot 2026-02-28 at 3.38.10 PM.png
alpha = 0.04;
beta = 1e4;
gamma = 3e7;
dydt = [ -alpha*y(1) + beta*y(2)*y(3);
alpha*y(1) – beta*y(2)*y(3) – gamma*y(2)^2;
gamma*y(2)^2 ];
end
function [t_all, u_all, error_history, h_history, stats] = …
dormand_prince_solver(f, tspan, u0, tol, h0)
% Initialization
t = tspan(1);
u = u0(:);
h = h0;
t_all = t;
u_all = u.’;
error_history = [];
h_history = [];
n_accepted = 0;
n_rejected = 0;
% —- Dormand–Prince 4(5) Coefficients —-
a = [1/5, 0, 0, 0, 0, 0;
3/40, 9/40, 0, 0, 0, 0;
44/45, -56/15, 32/9, 0, 0, 0;
19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0;
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0;
35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
b4 = [5179/57600, 0, 7571/16695, 393/640, …
-92097/339200, 187/2100, 1/40];
b5 = [35/384, 0, 500/1113, 125/192, …
-2187/6784, 11/84, 0];
c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
safety = 0.9;
h_min = 1e-12;
h_max = 1.0;
while t < tspan(2)
if h < h_min
error(‘Step size underflow: problem likely stiff.’);
end
if t + h > tspan(2)
h = tspan(2) – t;
end
% —- Stage calculations —-
k = zeros(length(u), 7);
k(:,1) = f(t, u);
for i = 2:7
k(:,i) = f(t + c(i)*h, …
u + h * (k(:,1:i-1) * a(i-1,1:i-1)’));
end
% —- 4th and 5th order solutions —-
y4 = u + h * (k * b4′);
y5 = u + h * (k * b5′);
error_val = norm(y5 – y4, inf);
if error_val <= tol
t = t + h;
u = y5;
t_all = [t_all; t];
u_all = [u_all; u.’];
error_history = [error_history; error_val];
h_history = [h_history; h];
n_accepted = n_accepted + 1;
else
n_rejected = n_rejected + 1;
end
if error_val == 0
h_new = 2*h;
else
h_new = h * safety * (tol / error_val)^(1/5);
end
h = min(max(h_new, h_min), h_max);
end
stats.accepted = n_accepted;
stats.rejected = n_rejected;
end I would like to know if this code can be improved or it well written? any comments will be greatly appreciated. Thanks
u0 = [1; 0; 0];
tspan = [0, 1e-3]; % small final time for explicit solver
tol = 1e-2;
h0 = 1e-6;
% Execute solver
[t, u, err, h_steps, stats] = dormand_prince_solver(@robertson_ode, tspan, u0, tol, h0);
% Part (a): Plot Solutions
figure;
subplot(2,1,1);
plot(t, u(:,1), ‘b’, t, u(:,3), ‘g’, ‘LineWidth’, 1.5);
legend(‘y1 (Reactant)’, ‘y3 (Product)’);
xlabel(‘Time t’); ylabel(‘Concentration’);
%title(‘Robertson Problem: Concentrations y1 and y3’);
grid on;
subplot(2,1,2);
semilogx(t, u(:,2), ‘r’, ‘LineWidth’, 1.5);
xlabel(‘Time t (Log Scale)’); ylabel(‘Concentration y2’);
%title(‘Robertson Problem: Intermediate Species y2’);
grid on;
% Part (b): Error and Stepsize Plots
figure;
subplot(2,1,1);
semilogy(t(2:end), err);
xlabel(‘Time t’); ylabel(‘||y_n^5 – y_n^4||_{infty} / h’);
%title(‘Part (b): Estimated Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps);
xlabel(‘Time t’); ylabel(‘Stepsize h’);
%title(‘Part (b): Stepsize h vs Time’);
grid on;
figure;
subplot(2,1,1);
semilogy(t(2:end), err, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Error’);
title(‘Estimated Local Error vs Time’);
grid on;
subplot(2,1,2);
semilogy(t(2:end), h_steps, ‘LineWidth’, 1.5);
xlabel(‘Time t’);
ylabel(‘Stepsize h’);
title(‘Stepsize vs Time’);
grid on;
function dydt = robertson_ode(~, y)
% Coefficients from Screenshot 2026-02-28 at 3.38.10 PM.png
alpha = 0.04;
beta = 1e4;
gamma = 3e7;
dydt = [ -alpha*y(1) + beta*y(2)*y(3);
alpha*y(1) – beta*y(2)*y(3) – gamma*y(2)^2;
gamma*y(2)^2 ];
end
function [t_all, u_all, error_history, h_history, stats] = …
dormand_prince_solver(f, tspan, u0, tol, h0)
% Initialization
t = tspan(1);
u = u0(:);
h = h0;
t_all = t;
u_all = u.’;
error_history = [];
h_history = [];
n_accepted = 0;
n_rejected = 0;
% —- Dormand–Prince 4(5) Coefficients —-
a = [1/5, 0, 0, 0, 0, 0;
3/40, 9/40, 0, 0, 0, 0;
44/45, -56/15, 32/9, 0, 0, 0;
19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0;
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0;
35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
b4 = [5179/57600, 0, 7571/16695, 393/640, …
-92097/339200, 187/2100, 1/40];
b5 = [35/384, 0, 500/1113, 125/192, …
-2187/6784, 11/84, 0];
c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1];
safety = 0.9;
h_min = 1e-12;
h_max = 1.0;
while t < tspan(2)
if h < h_min
error(‘Step size underflow: problem likely stiff.’);
end
if t + h > tspan(2)
h = tspan(2) – t;
end
% —- Stage calculations —-
k = zeros(length(u), 7);
k(:,1) = f(t, u);
for i = 2:7
k(:,i) = f(t + c(i)*h, …
u + h * (k(:,1:i-1) * a(i-1,1:i-1)’));
end
% —- 4th and 5th order solutions —-
y4 = u + h * (k * b4′);
y5 = u + h * (k * b5′);
error_val = norm(y5 – y4, inf);
if error_val <= tol
t = t + h;
u = y5;
t_all = [t_all; t];
u_all = [u_all; u.’];
error_history = [error_history; error_val];
h_history = [h_history; h];
n_accepted = n_accepted + 1;
else
n_rejected = n_rejected + 1;
end
if error_val == 0
h_new = 2*h;
else
h_new = h * safety * (tol / error_val)^(1/5);
end
h = min(max(h_new, h_min), h_max);
end
stats.accepted = n_accepted;
stats.rejected = n_rejected;
end code generation/improvement MATLAB Answers — New Questions
Graphs display upside down
I am using matlab online, 2021b
All graphs display upside down, including numbers, axis, and toolbox
Changing the renderer causes the graph to rerender and fixes the issue for graphs currently plotted, such as
set(gcf,’renderer’,’painters’)
set(gcf,’renderer’,’opengl’)
But does not stop new graphs from being upside down again.I am using matlab online, 2021b
All graphs display upside down, including numbers, axis, and toolbox
Changing the renderer causes the graph to rerender and fixes the issue for graphs currently plotted, such as
set(gcf,’renderer’,’painters’)
set(gcf,’renderer’,’opengl’)
But does not stop new graphs from being upside down again. I am using matlab online, 2021b
All graphs display upside down, including numbers, axis, and toolbox
Changing the renderer causes the graph to rerender and fixes the issue for graphs currently plotted, such as
set(gcf,’renderer’,’painters’)
set(gcf,’renderer’,’opengl’)
But does not stop new graphs from being upside down again. graph, graphics, gui, upside down MATLAB Answers — New Questions
How can I generate a License Report?
I am aware that I can use the license center to administer my licenses. I would like to see detailed information regarding license usage, is it possible to view this information in a comprehensive report?I am aware that I can use the license center to administer my licenses. I would like to see detailed information regarding license usage, is it possible to view this information in a comprehensive report? I am aware that I can use the license center to administer my licenses. I would like to see detailed information regarding license usage, is it possible to view this information in a comprehensive report? MATLAB Answers — New Questions
searchin for an M file in all folders of my PC.
I want to search for an M file in all the folders of my PC.
I do not want to search in every folder individually.I want to search for an M file in all the folders of my PC.
I do not want to search in every folder individually. I want to search for an M file in all the folders of my PC.
I do not want to search in every folder individually. file search MATLAB Answers — New Questions
Multiple curve fitting for parameter identification: multi-objective optimization
I have 4 plots of x (eps) and y (sigmaexp) and I would like to best-fit them in the equation given by sigmanum2. The criteria to best-fit is to minimize sum of squares i.e. least sum of squares defined by function fitness. I have done this successfully for 1 plot by using fminsearch i.e. 1 objective function. Since I have 4 plots, I have 4 objective functions (or 4 values in y i.e. [F1; F2; F3; F4]) which must be minimized simultaneously. How can I do this or which function must be used for multiobjective optimization? Please help. I have done research from my end in the relevant forums, however I am not able to do it. Here is my code.
s = table2array(ExpDataMatlab); %excel data I have imported. Attached for your reference.
eps = s(:,1); % represents x
sigmaexp = s(:,2); % represent y
epsdot1 = s(:,3);
A = 332.8668; %constant
B = 128.6184; %constant
n = 0.4234; %constant
%x = [2000, 10.2234];
%y = fitness(x,A,B,n,eps,epsdot1,sigmaexp); %I used this to check if my function defined at the
%bottom works properly or not
ObjFcn = @(x)fitness(x,A,B,n,eps,epsdot1,sigmaexp);
x0 = [12000 4];
options = optimset(‘Display’,’iter’,’PlotFcns’,@optimplotfval,’TolX’,1e-10,’TolFun’,1e-10,’Algorithm’,’interior-point’);
[sol, fval, exitflag, output] = fminsearch(ObjFcn,x0,options);
C = sol(1); P = sol(2);
sigmanum4 = (A+B*(eps.^n)).*(1+(epsdot1./C).^(1/P)); %deduced value after optimization
plot(s(1:99,1),s(1:99,2),’+m’);hold on;
plot(s(1:99,1),sigmanum4(1:99,1),’m’);
hold off;
function y = fitness(x,A,B,n,eps,epsdot1,sigmaexp)
sigmanum2 = (A+B*(eps.^n)).*(1+(epsdot1./x(1)).^(1/x(2)));
F1 = sum((sigmaexp(1:99,1)-sigmanum2(1:99,1)).^2);%values for 982 /s % represent residual sum of squares
%F2 = sum((sigmaexp(100:206,1)-sigmanum2(100:206,1)).^2);%values for 2000 /s
%F3 = sum((sigmaexp(207:296,1)-sigmanum2(207:296,1)).^2);%values for 2900 /s
%F4 = sum((sigmaexp(297:368,1)-sigmanum2(297:368,1)).^2);%values for 4000 /s
%y = [F1; F2; F3; F4]; % for multiobjective minimization
y = [F1]; %for single objective minimization
%y = (0.01042*F1+0.00962*F2+0.01111*F3+0.01389*F4); %weighted residual sum of squares (Milani 2008)
end
Thank you.I have 4 plots of x (eps) and y (sigmaexp) and I would like to best-fit them in the equation given by sigmanum2. The criteria to best-fit is to minimize sum of squares i.e. least sum of squares defined by function fitness. I have done this successfully for 1 plot by using fminsearch i.e. 1 objective function. Since I have 4 plots, I have 4 objective functions (or 4 values in y i.e. [F1; F2; F3; F4]) which must be minimized simultaneously. How can I do this or which function must be used for multiobjective optimization? Please help. I have done research from my end in the relevant forums, however I am not able to do it. Here is my code.
s = table2array(ExpDataMatlab); %excel data I have imported. Attached for your reference.
eps = s(:,1); % represents x
sigmaexp = s(:,2); % represent y
epsdot1 = s(:,3);
A = 332.8668; %constant
B = 128.6184; %constant
n = 0.4234; %constant
%x = [2000, 10.2234];
%y = fitness(x,A,B,n,eps,epsdot1,sigmaexp); %I used this to check if my function defined at the
%bottom works properly or not
ObjFcn = @(x)fitness(x,A,B,n,eps,epsdot1,sigmaexp);
x0 = [12000 4];
options = optimset(‘Display’,’iter’,’PlotFcns’,@optimplotfval,’TolX’,1e-10,’TolFun’,1e-10,’Algorithm’,’interior-point’);
[sol, fval, exitflag, output] = fminsearch(ObjFcn,x0,options);
C = sol(1); P = sol(2);
sigmanum4 = (A+B*(eps.^n)).*(1+(epsdot1./C).^(1/P)); %deduced value after optimization
plot(s(1:99,1),s(1:99,2),’+m’);hold on;
plot(s(1:99,1),sigmanum4(1:99,1),’m’);
hold off;
function y = fitness(x,A,B,n,eps,epsdot1,sigmaexp)
sigmanum2 = (A+B*(eps.^n)).*(1+(epsdot1./x(1)).^(1/x(2)));
F1 = sum((sigmaexp(1:99,1)-sigmanum2(1:99,1)).^2);%values for 982 /s % represent residual sum of squares
%F2 = sum((sigmaexp(100:206,1)-sigmanum2(100:206,1)).^2);%values for 2000 /s
%F3 = sum((sigmaexp(207:296,1)-sigmanum2(207:296,1)).^2);%values for 2900 /s
%F4 = sum((sigmaexp(297:368,1)-sigmanum2(297:368,1)).^2);%values for 4000 /s
%y = [F1; F2; F3; F4]; % for multiobjective minimization
y = [F1]; %for single objective minimization
%y = (0.01042*F1+0.00962*F2+0.01111*F3+0.01389*F4); %weighted residual sum of squares (Milani 2008)
end
Thank you. I have 4 plots of x (eps) and y (sigmaexp) and I would like to best-fit them in the equation given by sigmanum2. The criteria to best-fit is to minimize sum of squares i.e. least sum of squares defined by function fitness. I have done this successfully for 1 plot by using fminsearch i.e. 1 objective function. Since I have 4 plots, I have 4 objective functions (or 4 values in y i.e. [F1; F2; F3; F4]) which must be minimized simultaneously. How can I do this or which function must be used for multiobjective optimization? Please help. I have done research from my end in the relevant forums, however I am not able to do it. Here is my code.
s = table2array(ExpDataMatlab); %excel data I have imported. Attached for your reference.
eps = s(:,1); % represents x
sigmaexp = s(:,2); % represent y
epsdot1 = s(:,3);
A = 332.8668; %constant
B = 128.6184; %constant
n = 0.4234; %constant
%x = [2000, 10.2234];
%y = fitness(x,A,B,n,eps,epsdot1,sigmaexp); %I used this to check if my function defined at the
%bottom works properly or not
ObjFcn = @(x)fitness(x,A,B,n,eps,epsdot1,sigmaexp);
x0 = [12000 4];
options = optimset(‘Display’,’iter’,’PlotFcns’,@optimplotfval,’TolX’,1e-10,’TolFun’,1e-10,’Algorithm’,’interior-point’);
[sol, fval, exitflag, output] = fminsearch(ObjFcn,x0,options);
C = sol(1); P = sol(2);
sigmanum4 = (A+B*(eps.^n)).*(1+(epsdot1./C).^(1/P)); %deduced value after optimization
plot(s(1:99,1),s(1:99,2),’+m’);hold on;
plot(s(1:99,1),sigmanum4(1:99,1),’m’);
hold off;
function y = fitness(x,A,B,n,eps,epsdot1,sigmaexp)
sigmanum2 = (A+B*(eps.^n)).*(1+(epsdot1./x(1)).^(1/x(2)));
F1 = sum((sigmaexp(1:99,1)-sigmanum2(1:99,1)).^2);%values for 982 /s % represent residual sum of squares
%F2 = sum((sigmaexp(100:206,1)-sigmanum2(100:206,1)).^2);%values for 2000 /s
%F3 = sum((sigmaexp(207:296,1)-sigmanum2(207:296,1)).^2);%values for 2900 /s
%F4 = sum((sigmaexp(297:368,1)-sigmanum2(297:368,1)).^2);%values for 4000 /s
%y = [F1; F2; F3; F4]; % for multiobjective minimization
y = [F1]; %for single objective minimization
%y = (0.01042*F1+0.00962*F2+0.01111*F3+0.01389*F4); %weighted residual sum of squares (Milani 2008)
end
Thank you. multiple curve fitting, parameter identification, multi-objective optimization MATLAB Answers — New Questions
Maltab Version R2015a required for student download
Hello,
I am a student from Heidelberg University. For my master thesis I would need the Matlab Version R2015a or older.
I could not find this version in the student download version (Free version).
Thank you for your help.
Kind regards,
Max BickelhauptHello,
I am a student from Heidelberg University. For my master thesis I would need the Matlab Version R2015a or older.
I could not find this version in the student download version (Free version).
Thank you for your help.
Kind regards,
Max Bickelhaupt Hello,
I am a student from Heidelberg University. For my master thesis I would need the Matlab Version R2015a or older.
I could not find this version in the student download version (Free version).
Thank you for your help.
Kind regards,
Max Bickelhaupt matlab version r 2015a MATLAB Answers — New Questions
How can I store data from my simulink model in and and use it during runtime.
How can I store data from my simulink model in and and use it during runtime. I used the data store blocks and trigger subsystem to enable data store only during specific events. But at a sample time of 1e-4 sec there’s a lot of data. I want to only buffer 50 points.How can I store data from my simulink model in and and use it during runtime. I used the data store blocks and trigger subsystem to enable data store only during specific events. But at a sample time of 1e-4 sec there’s a lot of data. I want to only buffer 50 points. How can I store data from my simulink model in and and use it during runtime. I used the data store blocks and trigger subsystem to enable data store only during specific events. But at a sample time of 1e-4 sec there’s a lot of data. I want to only buffer 50 points. data store simulink MATLAB Answers — New Questions
dot product for complex vector
Hello,
In the Matlab example, you have the dot product of the following two vectors A and B and its answer is vector C.
A = [1+i 1-i -1+i -1-i];
B = [3-4i 6-2i 1+2i 4+3i];
Calculate the dot product of A and B.
C = dot(A,B)
C = 1.0000 – 5.0000i
However, when I calculate it, I have vector C = 7 – 17i
That is, I have C vector results as follows below
(1+i) * (3-4i) + (1-i) * (6-2i) + (-1+i) * (1+2i) + (-1-i) * (4+3i) =
(7-i) +( 4-8i) + (-3-i) + (-1-7i) =
7 – 17i.
Hence, could you please tell me how the Matlab got the results (or show me manually how Matlab got the dot product answer) as I have different results than Matlab calculated using dot product?
Thank you,
CharlesHello,
In the Matlab example, you have the dot product of the following two vectors A and B and its answer is vector C.
A = [1+i 1-i -1+i -1-i];
B = [3-4i 6-2i 1+2i 4+3i];
Calculate the dot product of A and B.
C = dot(A,B)
C = 1.0000 – 5.0000i
However, when I calculate it, I have vector C = 7 – 17i
That is, I have C vector results as follows below
(1+i) * (3-4i) + (1-i) * (6-2i) + (-1+i) * (1+2i) + (-1-i) * (4+3i) =
(7-i) +( 4-8i) + (-3-i) + (-1-7i) =
7 – 17i.
Hence, could you please tell me how the Matlab got the results (or show me manually how Matlab got the dot product answer) as I have different results than Matlab calculated using dot product?
Thank you,
Charles Hello,
In the Matlab example, you have the dot product of the following two vectors A and B and its answer is vector C.
A = [1+i 1-i -1+i -1-i];
B = [3-4i 6-2i 1+2i 4+3i];
Calculate the dot product of A and B.
C = dot(A,B)
C = 1.0000 – 5.0000i
However, when I calculate it, I have vector C = 7 – 17i
That is, I have C vector results as follows below
(1+i) * (3-4i) + (1-i) * (6-2i) + (-1+i) * (1+2i) + (-1-i) * (4+3i) =
(7-i) +( 4-8i) + (-3-i) + (-1-7i) =
7 – 17i.
Hence, could you please tell me how the Matlab got the results (or show me manually how Matlab got the dot product answer) as I have different results than Matlab calculated using dot product?
Thank you,
Charles dot product MATLAB Answers — New Questions
Microsoft and OpenAI joint statement on continuing partnership
Since 2019, Microsoft and OpenAI have worked together to advance artificial intelligence responsibly and make its benefits broadly accessible. What began as a research partnership has grown into one of the most consequential collaborations in technology — grounded in mutual trust, deep technical integration, and a long‑term commitment to innovation.
As conversations around AI investments and partnerships grow and as OAI announces new funding and new partners as they did today, we want to ensure these announcements are understood within the existing construct of our partnership. Nothing about today’s announcements in any way changes the terms of the Microsoft and OpenAI relationship that have been previously shared in our joint blog in October 2025.
The partnership remains strong and central. Microsoft and OpenAI continue to work closely across research, engineering, and product development, building on years of deep collaboration and shared success.
Our IP relationship continues unchanged. Microsoft maintains its exclusive license and access to intellectual property across OpenAI models and products. Collaborations like the partnership between OpenAI and Amazon were always contemplated under our agreements and Microsoft is excited to see what they build together.
Our commercial and revenue share relationship remains unchanged. The ongoing revenue share arrangement remains unchanged and has always included sharing revenue from partnerships between OpenAI and other cloud providers.
Azure remains the exclusive cloud provider of stateless OpenAI APIs. Microsoft is the exclusive cloud provider for stateless APIs that provide access to OpenAI’s models and IP. These APIs can be purchased from Microsoft or directly from OpenAI. Customers and developers benefit from Azure’s global infrastructure, security, and enterprise-grade capabilities at scale. Any stateless API calls to OpenAI models that result from a collaboration between OpenAI and any third party – including Amazon – would be hosted on Azure.
OpenAI’s first party products, including Frontier, will continue to be hosted on Azure.
AGI definition and processes are unchanged. The contractual definition of AGI and the process for determining if it has been achieved remains the same.
The partnership supports OpenAI’s growth. As OpenAI scales, it continues to have flexibility to commit to additional compute elsewhere, including through large-scale infrastructure initiatives such as the Stargate project.
The partnership was designed to give Microsoft and OpenAI room to pursue new opportunities independently, while continuing to collaborate, which each company is doing, together and independently.
We remain committed to our partnership and to the shared mission that brought us together. We continue to work side‑by‑side to deliver powerful AI tools, advance responsible development, and ensure that AI benefits people and organizations everywhere.
The post Microsoft and OpenAI joint statement on continuing partnership appeared first on The Official Microsoft Blog.
Since 2019, Microsoft and OpenAI have worked together to advance artificial intelligence responsibly and make its benefits broadly accessible. What began as a research partnership has grown into one of the most consequential collaborations in technology — grounded in mutual trust, deep technical integration, and a long‑term commitment to innovation. As conversations around AI investments and partnerships…
The post Microsoft and OpenAI joint statement on continuing partnership appeared first on The Official Microsoft Blog.Read More
Have installed Matlab 2025b on my Macbook air M2 but it crashes at startup. Have tried a new installation but didn’t work.
Hello there, I need help installing matlab 2025b on my macbook running Tahoe 26.1. The instalation runs smothly but once instaled, Matlab crashes at startup. I have erased and tried a fresh instalation but have not worked. What would you suggest me to do?
Yours sincerely,
DanielHello there, I need help installing matlab 2025b on my macbook running Tahoe 26.1. The instalation runs smothly but once instaled, Matlab crashes at startup. I have erased and tried a fresh instalation but have not worked. What would you suggest me to do?
Yours sincerely,
Daniel Hello there, I need help installing matlab 2025b on my macbook running Tahoe 26.1. The instalation runs smothly but once instaled, Matlab crashes at startup. I have erased and tried a fresh instalation but have not worked. What would you suggest me to do?
Yours sincerely,
Daniel matlab installation, mac, 2025b, tahoe, crash MATLAB Answers — New Questions
Please help me to solve this simple error in definition of dy(3)
proj()
function sol= proj
myLegend1 = {};
myLegend2 = {};
k0=386; ce=3.831*10^2;rho=89.54*10^2;alfat=1.78*10^-5;taw=0.5;Tnot=2.93*10^2;mu=38.6*10^5; lamda=77.6*10^5;
%lamda=77.6e^9;
%mu=38.6*e^9;lamda=77.6e+9
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0.3 ];
for i =1:numel(rr)
a= rr(i);
s=0.01;h=0.01;b=0.01;
y0 = [0,0,0,0,0,0,1,0,0];
disp(a7)
options =bvpset(‘stats’,’on’,’RelTol’,1e-4);
m = linspace(0,2);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(sol.x,abs(sol.y(1,:)))
grid on,hold on
myLegend1{i}=[‘alfa= ‘,num2str(rr(i))];
figure(2)
plot(sol.x,(sol.y(7,:)))
title(‘Temperature’)
grid on,hold on
myLegend2{i}=[‘alfa= ‘,num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
function dy = projfun(x,y)
dy = zeros(9,1);
E = y(1);
dE = y(2);
ddE=y(3);
u = y(4);
du = y(5);
ddu=y(6);
t = y(7);
dt = y(8);
ddt=y(9);
dy(1) = dE;
dy(2)=(1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du));
dy(3)=((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)^2))*(((((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)))*((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*s*ddt-a3*dt-a*(a4*(s+h)+a1*(s+h))*dt*dE-a*(s^2+a1*h^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+dt*dE)))+(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(s^2+a1*h^2-(x+b+1))*(1-a*t)*((2*a*a3*h)*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*dt-ddt*a3*h-a*(a4+a1)*(s+h)*dt*du-a*(h^2+a1*s^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+ddt*dE)-a*2*(h+s)*dt*dE-a*s*h*(a4+a1)*(dt*((1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-1*a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))))+dt*du)-(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)*(2*s+2*a1*h-a*(s^2+a1*h^2-(x+b+1)^2)*dt-(2*a*(s+a1*h))*t))-((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*(a*s*h*(a4+a1))*dt*dE)-(s^2+a1*h^2-(x+b+1)^2)*(1-a*t)*(2*a*a3*h*t*dt-a3*h*dt-a*(s^2+a1*h^2)*dt*dE-a*s*h*(a4+a1)*dt*du))*((s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2)+2*a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2)*t*dt+t^2*(a*(s^2+a1*h^2-(x+b+1)^2)*a*2*(h+a1*s)+a*2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2))-t*((s^2+a1*h^2-(x+b+1)^2)*a*2*(h+s*a1)+2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2)+a*(s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+a*2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2))-2*a2*s*h(1-a*t)*(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)-dt*(a*(s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))));
dy(4)=du;
dy(5)=(1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du))));
dy(6) =(1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE));
dy(7)=dt;
dy(8)=ddt;
dy(9)=(1/((a*(x+b+1)*(a5*(s^2+h^2)-a8*(x+b+1)*(x+b+1)))*t-((x+b+1)*(a5*s^2+a5*h^2-a8*(x+b+1)*(x+b+1)))))*((s^2+h^2+2*a5*s*+2*a5*h-a6*(x+b+1)*(x+b+1))*ddt+((2*a*a7*(x+b+1))*t-2*(x+b+1)*(a7-a*a11))*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+((a*2*a7)*t-(2*(x+b+1)*(a7-2*a*a11)))*((1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)))-a*(s^2+h^2+2*a5*s+2*a5*h-a6*(x+b+1)*(x+b+1))*t*ddt-a*(s^2+h^2+a5*s+a5*h)*dt*dt-a*(x+b+1)*(a5*s^2+a5*h^2+a10*(x+b+1)*(x+b+1))*dt*ddt);
end
function res = projbc(ya,yb)
res = [ya(1);
ya(4);
ya(7)-1;
ya(3);
ya(5);
yb(1);
yb(4);
yb(7)-0.5;
yb(9)];
end
endproj()
function sol= proj
myLegend1 = {};
myLegend2 = {};
k0=386; ce=3.831*10^2;rho=89.54*10^2;alfat=1.78*10^-5;taw=0.5;Tnot=2.93*10^2;mu=38.6*10^5; lamda=77.6*10^5;
%lamda=77.6e^9;
%mu=38.6*e^9;lamda=77.6e+9
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0.3 ];
for i =1:numel(rr)
a= rr(i);
s=0.01;h=0.01;b=0.01;
y0 = [0,0,0,0,0,0,1,0,0];
disp(a7)
options =bvpset(‘stats’,’on’,’RelTol’,1e-4);
m = linspace(0,2);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(sol.x,abs(sol.y(1,:)))
grid on,hold on
myLegend1{i}=[‘alfa= ‘,num2str(rr(i))];
figure(2)
plot(sol.x,(sol.y(7,:)))
title(‘Temperature’)
grid on,hold on
myLegend2{i}=[‘alfa= ‘,num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
function dy = projfun(x,y)
dy = zeros(9,1);
E = y(1);
dE = y(2);
ddE=y(3);
u = y(4);
du = y(5);
ddu=y(6);
t = y(7);
dt = y(8);
ddt=y(9);
dy(1) = dE;
dy(2)=(1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du));
dy(3)=((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)^2))*(((((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)))*((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*s*ddt-a3*dt-a*(a4*(s+h)+a1*(s+h))*dt*dE-a*(s^2+a1*h^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+dt*dE)))+(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(s^2+a1*h^2-(x+b+1))*(1-a*t)*((2*a*a3*h)*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*dt-ddt*a3*h-a*(a4+a1)*(s+h)*dt*du-a*(h^2+a1*s^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+ddt*dE)-a*2*(h+s)*dt*dE-a*s*h*(a4+a1)*(dt*((1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-1*a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))))+dt*du)-(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)*(2*s+2*a1*h-a*(s^2+a1*h^2-(x+b+1)^2)*dt-(2*a*(s+a1*h))*t))-((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*(a*s*h*(a4+a1))*dt*dE)-(s^2+a1*h^2-(x+b+1)^2)*(1-a*t)*(2*a*a3*h*t*dt-a3*h*dt-a*(s^2+a1*h^2)*dt*dE-a*s*h*(a4+a1)*dt*du))*((s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2)+2*a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2)*t*dt+t^2*(a*(s^2+a1*h^2-(x+b+1)^2)*a*2*(h+a1*s)+a*2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2))-t*((s^2+a1*h^2-(x+b+1)^2)*a*2*(h+s*a1)+2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2)+a*(s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+a*2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2))-2*a2*s*h(1-a*t)*(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)-dt*(a*(s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))));
dy(4)=du;
dy(5)=(1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du))));
dy(6) =(1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE));
dy(7)=dt;
dy(8)=ddt;
dy(9)=(1/((a*(x+b+1)*(a5*(s^2+h^2)-a8*(x+b+1)*(x+b+1)))*t-((x+b+1)*(a5*s^2+a5*h^2-a8*(x+b+1)*(x+b+1)))))*((s^2+h^2+2*a5*s*+2*a5*h-a6*(x+b+1)*(x+b+1))*ddt+((2*a*a7*(x+b+1))*t-2*(x+b+1)*(a7-a*a11))*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+((a*2*a7)*t-(2*(x+b+1)*(a7-2*a*a11)))*((1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)))-a*(s^2+h^2+2*a5*s+2*a5*h-a6*(x+b+1)*(x+b+1))*t*ddt-a*(s^2+h^2+a5*s+a5*h)*dt*dt-a*(x+b+1)*(a5*s^2+a5*h^2+a10*(x+b+1)*(x+b+1))*dt*ddt);
end
function res = projbc(ya,yb)
res = [ya(1);
ya(4);
ya(7)-1;
ya(3);
ya(5);
yb(1);
yb(4);
yb(7)-0.5;
yb(9)];
end
end proj()
function sol= proj
myLegend1 = {};
myLegend2 = {};
k0=386; ce=3.831*10^2;rho=89.54*10^2;alfat=1.78*10^-5;taw=0.5;Tnot=2.93*10^2;mu=38.6*10^5; lamda=77.6*10^5;
%lamda=77.6e^9;
%mu=38.6*e^9;lamda=77.6e+9
c0=sqrt((lamda+2*mu)/(rho)); Betanot=(3*lamda+2*mu)*alfat; a1=mu/(lamda+2*mu);a2=(mu+lamda)/(lamda+2*mu);a3=(Betanot*Tnot)/(lamda+2*mu);omega=(rho*ce)/(k0);
a4=lamda/(lamda+2*mu);a5=(k0*omega*c0^2)/(k0);a6=(rho*ce*c0^2)/(k0);
a7=(Betanot*c0^2)/(k0); a8=a6*taw; a9=a7*taw; a10=rho*ce*taw*omega*c0^4/(k0); a11=Betanot*taw*omega*c0^4/(k0);w=rho*ce/(k0);
rr = [0.3 ];
for i =1:numel(rr)
a= rr(i);
s=0.01;h=0.01;b=0.01;
y0 = [0,0,0,0,0,0,1,0,0];
disp(a7)
options =bvpset(‘stats’,’on’,’RelTol’,1e-4);
m = linspace(0,2);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
figure(1)
plot(sol.x,abs(sol.y(1,:)))
grid on,hold on
myLegend1{i}=[‘alfa= ‘,num2str(rr(i))];
figure(2)
plot(sol.x,(sol.y(7,:)))
title(‘Temperature’)
grid on,hold on
myLegend2{i}=[‘alfa= ‘,num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
figure(2)
legend(myLegend2)
function dy = projfun(x,y)
dy = zeros(9,1);
E = y(1);
dE = y(2);
ddE=y(3);
u = y(4);
du = y(5);
ddu=y(6);
t = y(7);
dt = y(8);
ddt=y(9);
dy(1) = dE;
dy(2)=(1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du));
dy(3)=((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)^2))*(((((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2)))*((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*s*ddt-a3*dt-a*(a4*(s+h)+a1*(s+h))*dt*dE-a*(s^2+a1*h^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+dt*dE)))+(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(s^2+a1*h^2-(x+b+1))*(1-a*t)*((2*a*a3*h)*(t*ddt+dt*dt)+2*a*a3*t*dt-a3*dt-ddt*a3*h-a*(a4+a1)*(s+h)*dt*du-a*(h^2+a1*s^2)*(dt*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))+ddt*dE)-a*2*(h+s)*dt*dE-a*s*h*(a4+a1)*(dt*((1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-1*a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du)))))+dt*du)-(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)*(2*s+2*a1*h-a*(s^2+a1*h^2-(x+b+1)^2)*dt-(2*a*(s+a1*h))*t))-((a2*s*h-a*a2*s*h*t)*(2*a*a3*s*t*dt-a3*s*dt-a*(s^2+a1*h^2)*dt*du-a*(a*s*h*(a4+a1))*dt*dE)-(s^2+a1*h^2-(x+b+1)^2)*(1-a*t)*(2*a*a3*h*t*dt-a3*h*dt-a*(s^2+a1*h^2)*dt*dE-a*s*h*(a4+a1)*dt*du))*((s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2)+2*a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2)*t*dt+t^2*(a*(s^2+a1*h^2-(x+b+1)^2)*a*2*(h+a1*s)+a*2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2))-t*((s^2+a1*h^2-(x+b+1)^2)*a*2*(h+s*a1)+2*(s+a1*h)*a*(h^2+a1*s^2-(x+b+1)^2)+a*(s^2+a1*h^2-(x+b+1)^2)*2*(h+a1*s)+a*2*(s+a1*h)*(h^2+a1*s^2-(x+b+1)^2))-2*a2*s*h(1-a*t)*(a2*(s+h)-a*a2*s*h*dt-a*a2*(s+h)*t)-dt*(a*(s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))));
dy(4)=du;
dy(5)=(1/((a*(s^2+a1*h^2-(x+b+1)^2)*t-(s^2+a1*h^2-(x+b+1)^2))))*((-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)+(a2*s*h-a*a2*s*h*t)*((1/((s^2+a1*h^2-(x+b+1)^2)*(h^2+a1*s^2-(x+b+1)^2)+(a*(s^2+a1*h^2-(x+b+1)^2)*a*(h^2+a1*s^2-(x+b+1)^2))*t^2-(a*(s^2+a1*h^2-(x+b+1)^2)*((s^2+a1*h^2-(x+b+1)^2))+(s^2+a1*h^2-(x+b+1)^2)*a*(s^2+a1*h^2-(x+b+1)^2))*t-(a2*s*h-a*a2*s*h*t)^2))*((a2*s*h-a*a2*s*h*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-(a*a4*s*h+a*a1*s*h)*dt*dE)-((s^2+a1*h^2-(x+b+1)^2)-a*((s^2+a1*h^2-(x+b+1)^2))*t)*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*(a4*s*h+a1*s*h)*dt*du))));
dy(6) =(1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1).^2))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE));
dy(7)=dt;
dy(8)=ddt;
dy(9)=(1/((a*(x+b+1)*(a5*(s^2+h^2)-a8*(x+b+1)*(x+b+1)))*t-((x+b+1)*(a5*s^2+a5*h^2-a8*(x+b+1)*(x+b+1)))))*((s^2+h^2+2*a5*s*+2*a5*h-a6*(x+b+1)*(x+b+1))*ddt+((2*a*a7*(x+b+1))*t-2*(x+b+1)*(a7-a*a11))*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+((a*2*a7)*t-(2*(x+b+1)*(a7-2*a*a11)))*((1/((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*((a2*s*h-(a*a2*s*h)*t)*((1/((((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t)*(((s^2+a1*h^2-(x+b+1)*(x+b+1)))-a*(s^2+a1*h^2-(x+b+1)*(x+b+1))*t))-(a2*s*h-(a*a2*s*h)*t)^2))*((a2*s*h-(a*a2*s*h)*t)*(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)-(((s^2+a1*h^2-(x+b+1)*(x+b+1))-(a*(s^2+a1*h^2-(x+b+1)*(x+b+1)))*t))*(-a3*h*dt+2*a*a3*h*t*dt-a*(h^2+a1*s^2)*dt*dE-a*s*h*(a4+a1)*dt*du)))+(-a3*s*dt+2*a*a3*s*t*dt-a*(s^2+a1*h^2)*dt*du-a*s*h*(a4+a1)*dt*dE)))-a*(s^2+h^2+2*a5*s+2*a5*h-a6*(x+b+1)*(x+b+1))*t*ddt-a*(s^2+h^2+a5*s+a5*h)*dt*dt-a*(x+b+1)*(a5*s^2+a5*h^2+a10*(x+b+1)*(x+b+1))*dt*ddt);
end
function res = projbc(ya,yb)
res = [ya(1);
ya(4);
ya(7)-1;
ya(3);
ya(5);
yb(1);
yb(4);
yb(7)-0.5;
yb(9)];
end
end pvb4c MATLAB Answers — New Questions
Simulation issues with current control with a 3 phase inverter
Im working on the following simulation.
Where im using the programmable voltage source as a representation of the grid and than im using a inductance and resistance with 3mH and 50 mOhm and than im using a Average value voltage source converter representing a bidirectional inverter and than a dc link with a batterie with 800V.
The following is my current control algorithm,
The control signal is than put into a PWM Timing and Waveform Generator to control the inverter.
The issue im currently facing is that i cant control the current for some reason the current id is not following the id* current, i tried tuning the PI controller but that doesnt help. I still get the following currents even tho they should go to 0 since id* and iq* are both 0.
I then tried to run the vd voltage as the vd inv voltage and the vq as the vq inv Voltage but there is still some residual current in the system even tho it should be 0 at that point.
I really tried a lot of different approaches like measuring the current between inductance and inverter, and also played with the sample time but i never get the current id and iq to follow id* and iq* i always get this current form and i dont really understand why.
Is my simulation wrong physically or am i not seeing something?
I also put a copy of my simulation below in case someone wants to test with it.
Any kind of help or new leads as why this current controller doesnt work would be much appreciated.Im working on the following simulation.
Where im using the programmable voltage source as a representation of the grid and than im using a inductance and resistance with 3mH and 50 mOhm and than im using a Average value voltage source converter representing a bidirectional inverter and than a dc link with a batterie with 800V.
The following is my current control algorithm,
The control signal is than put into a PWM Timing and Waveform Generator to control the inverter.
The issue im currently facing is that i cant control the current for some reason the current id is not following the id* current, i tried tuning the PI controller but that doesnt help. I still get the following currents even tho they should go to 0 since id* and iq* are both 0.
I then tried to run the vd voltage as the vd inv voltage and the vq as the vq inv Voltage but there is still some residual current in the system even tho it should be 0 at that point.
I really tried a lot of different approaches like measuring the current between inductance and inverter, and also played with the sample time but i never get the current id and iq to follow id* and iq* i always get this current form and i dont really understand why.
Is my simulation wrong physically or am i not seeing something?
I also put a copy of my simulation below in case someone wants to test with it.
Any kind of help or new leads as why this current controller doesnt work would be much appreciated. Im working on the following simulation.
Where im using the programmable voltage source as a representation of the grid and than im using a inductance and resistance with 3mH and 50 mOhm and than im using a Average value voltage source converter representing a bidirectional inverter and than a dc link with a batterie with 800V.
The following is my current control algorithm,
The control signal is than put into a PWM Timing and Waveform Generator to control the inverter.
The issue im currently facing is that i cant control the current for some reason the current id is not following the id* current, i tried tuning the PI controller but that doesnt help. I still get the following currents even tho they should go to 0 since id* and iq* are both 0.
I then tried to run the vd voltage as the vd inv voltage and the vq as the vq inv Voltage but there is still some residual current in the system even tho it should be 0 at that point.
I really tried a lot of different approaches like measuring the current between inductance and inverter, and also played with the sample time but i never get the current id and iq to follow id* and iq* i always get this current form and i dont really understand why.
Is my simulation wrong physically or am i not seeing something?
I also put a copy of my simulation below in case someone wants to test with it.
Any kind of help or new leads as why this current controller doesnt work would be much appreciated. current control, error, inverter control, bidirectional grid simulation, simulation MATLAB Answers — New Questions
